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1.  INTRODUCTION 

In  circuit  theory,  the  impulse  response  that  characterizes  a linear 
circuit  may  be  determined  from  the  knowledge  of  the  singularities  of  the 
response  function  in  the  complex  frequency  plane  and  the  corresponding 
residues.  The  impulse  response  of  the  circuit  is  then  simply  a summation 
of  all  the  residues  multiplied  by  exponentially  damped  sinusoids.  In 
recent  years  a similar  approach  has  been  used  to  characterize  the  tran- 
sient phenomenon  of  electromagnetic  scattering  and  antenna  problems.  This 
approach  is  known  as  the  singularity  expansion  method  (SEM)  [1],  [2]. 

The  development  of  the  singularity  expansion  method  arose  from  in- 
sight into  the  general  characteristics  of  typical  transient  response 
behavior  observed  on  various  electromagnetic  structures.  The  transient 
response  waveforms  of  these  structures  appear  to  be  dominated  by  a few 
exponentially  damped  sinusoids.  This  observation  is  especially  apparent 
if  one  looks  at  the  late  time  response  of  slender  or  thin  wire  structures. 
The  resonant  frequency,  damping  constant,  and  current  distribution  of 
some  of  these  resonant  modes  had  been  calculated  for  the  thin  wire  and 
prolate  spheroid  as  early  as  1930  [3],  [4].  However,  not  until  the  sin- 
gularity expansion  method  came  into  being  was  it  possible  to  determine  the 
modal  resonances  and  the  excitation  coefficient  of  each  mode  for  a 
structure  with  an  arbitrary  incident  wave. 

In  the  singularity  expansion  method  the  transient  response  is  written 
as  a sum  of  exponentially  damped  sinusoids.  In  order  to  write  this  sum, 
it  is  necessary  to  first  determine  the  location  of  the  complex  natural 
frequencies  or  pole  singularities  of  the  structure  being  studied.  The 
conventional  approach  for  determining  the  singularities  of  a system  is 
based  on  an  iterative  search  procedure  that  seeks  the  zeros  of  the  system 
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determinant  in  the  complex  frequency  plane.  This  approach  has  been  used 
successfully  by  many  people  and  has  given  extremely  satisfying  results. 

A few  illustrations  are  given  in  References  [5]  - [7].  Since  the  iter- 
ative search  is  a slow  procedure,  it  is  usually  economical  to  extract 
only  a few  poles.  Also  in  the  conventional  techniques,  the  poles  cannot 
be  extracted  from  the  time-domain  formulation  of  the  problem. 

Recent  development  of  methods  for  the  direct  production  of  both 
numerical  [8],  [9]  and  experimental  [10],  [11]  transient  electromagnetic 
response  data  has  generated  considerable  interest  in  the  possibility  of 
direct  extraction  of  the  poles  and  their  accompanying  residues  from 
given  time-domain  system  response.  This  thesis  presents  a numerical 
technique  for  systematically  determining  the  system  singularities  from 
the  transient  response  data  of  that  system.  The  approach  that  will  be 
developed  here  is  based  on  Prony's  algorithm  which  was  first  published 
in  1795  [12]  and  has  appeared  in  a few  good  numerical  analysis  texts 
[13],  [14]. 

In  Chapter  2 the  mathematical  notation  is  developed  for  the  impulse 
response  function  for  electromagnetic  antennas  and  scatterers.  The 
assumptions  used  in  reducing  the  impulse  response  to  a sum  of  exponentials 
are  also  presented.  The  necessity  of  removing  the  influence  of  the 
driving  function  from  the  transient  response  if  the  driving  function  is 
not  a true  delta  function  is  discussed  in  detail. 

Prony's  method  is  developed  in  detail  in  Chapter  3.  The  method  is 
applicable  to  systems  containing  multiple  as  well  as  simple  poles.  Sev- 
eral numerical  examples  are  presented  and  the  results  are  analyzed. 

These  results  are  used  to  establish  guidelines  for  the  use  of  the  method 
and  to  bring  out  some  of  the  problems  associated  with  the  method. 
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A systematic  procedure  by  which  the  number  of  poles  inherent  in  the 
transient  response  can  be  determined  is  necessary  in  order  to  properly 
use  Prony's  method.  Chapter  4 presents  two  techniques  for  doing  this. 

The  first  technique  is  based  on  an  orthogonalization  procedure  and  the 
second  technique  is  based  on  an  eigenvalue  approach.  Both  methods  have 
advantages  and  disadvantages  inherent  in  their  implementation.  These 
advantages  and  disadvantages  are  presented  by  showing  several  numerical 
examples. 

Tt  is  shown  in  Chapter  5 that  noise  seriously  affects  the  poles 
returned  by  Prony's  method.  In  some  cases  the  noise  level  is  high  enough 
to  completely  corrupt  the  results.  Several  statistical  studies  are  pre- 
sented which  relate  the  standard  deviation  of  the  noise  to  the  quality  of 
the  results  obtained  from  Prony's  method. 

Chapter  6 discusses  several  of  the  applications  in  which  the  de- 
veloped methods  can  be  implemented.  Some  of  these  topics  are:  system 

analysis,  radar  target  recognition,  the  study  of  spectral  characteristics, 
and  data  reduction  and  extrapolation.  Examples  are  used  to  show  that 
Prony's  method  applied  to  these  problems  has  definite  numerical  advan- 
tages over  the  conventional  approaches  used. 

An  alternative  to  Prony's  method  is  presented  in  Chapter  7.  This 
approach  is  the  Pade  approximation  but  it  is  extremely  limited  in  its 
usefulness.  Chapter  8 presents  conclusions  and  recommendations  for 


further  study. 


2.  MATHEMATICAL  FORM  OF  THE  TRANSIENT  RESPONSE 

This  chapter  introduces  the  mathematical  form  of  the  impulse  re- 
sponse function  for  an  electromagnetic  scatterer  or  antenna.  The  impulse 
response  will  be  written  in  both  the  s-plane  (Laplace  transform)  domain 
and  in  the  time  domain.  The  assumptions  necessary  for  reducing  the  impulse 
response  to  a simple  sum  of  exponentials  are  discussed.  It  is  necessary 
to  reduce  the  response  to  a sum  of  exponentials  so  that  the  methods  of 
Chapter  3 may  be  applied.  The  form  of  the  response  functions  for  an 
arbitrary  excitation  is  also  studied. 

2 . 1 The  Impulse  Response 

The  normalized  s-plane  impulse  or  delta  function  response  for  a 
finite-sized  object  that  has  only  pole  singularities  in  free  space  is 
generally  written  [11,  [2] 

“ -m. 

H(r,s)  = l n.(s  ,p)  v (r)  (s  - s ) 1 +W(r,s,p)  (2.1) 

i=l 

where  the  above  terms  are  defined  as: 

s^  - natural  frequency,  pole  singularity,  natural  resonances. 

This  is  a complex  frequency  for  which  the  system  has  a re- 
sponse when  no  forcing  function  is  applied.  The  poles  must 
appear,  in  complex  conjugate  pairs  or  lie  on  the  real  axis. 

The  poles  also  must  lie  in  the  negative  half  of  the  s-plane. 
v^(r)  - natural  mode. 

This  is  the  response  of  the  system  at  s^  which  depends  on 


the  position  r on  the  structure  and  the  object  parameters 


r\ i (s £ » P ) ~ coupling  coefficient. 

This  is  the  strength  of  the  natural  oscillation  s^  in  terms 
of  the  system  and  the  incident  wave  parameters.  It  is  inde- 
pendent of  position. 

p - polarization  of  incident  plane  wave, 

r - the  position  vector. 

» 

This  is  the  position  on  the  structure  at  which  the  transient 
response  is  being  measured  or  observed. 

t*  Vi 

nu  - the  multiplicity  of  the  i pole. 

The  term  W(r,s,p)  is  an  entire  function  of  s and  dependent  on  the  form 
of  the  coupling  coefficient  and  the  incident  wave.  In  the  most  general 
case  this  term  is  required  by  the  Mittag-Lef f ler  Theorem  [15]  in  order 
to  guarantee  convergence  of  the  infinite  series.  It  has  been  hypothe- 
sized that  the  entire  function  is  not  needed  in  most  electromagnetics 
problems  and  could  normally  be  neglected. 


The  impulse  response  (2.1)  can  be  written  in  the  time  domain  as 
* s t 

H(r,t)  = u ( t - t„)  l n.(s  ,p)  v '“)  e 1 (2.2) 

i=l  1 1 

where  the  entire  function  has  been  neglected  and  all  poles  have  been 
assumed  to  be  simple.  The  step  function  u(t  - tg)  is  present  so  that 
the  response  does  not  start  until  tg,  the  time  at  which  the  response 
begins  at  the  particular  observation  position  r on  the  body.  Since 
the  entire  function  has  been  neglected,  it  is  necessary  to  require  that 
tg  be  greater  than  zero  because  a delta  function  source  applied  at  r 
at  t * 0 yields  a delta  function  at  t * 0 in  the  impulse  response  which 
cannot  be  represented  by  the  exponential  terms. 


As  an  example,  consider  a perfectly  conducting  finite  dipole  driven 


at  t » 0 by  a delta  function  plane  wave  with  the  E field  polarized  in 
the  direction  of  the  dipole  axis.  The  impulse  response  of  the  induced 
current  at  any  position  r on  the  dipole  is  of  the  form 

s.t 


I(r,t) 


l A (r)  e 1 + B . (r)  5(t)  u(t)  . 

i=l  1 


(2.3) 


If  the  induced  current  is  expressed  as  a sum  of  complex  exponentials  it 
is  written  as 


Kr.t) 


' 00  _ s.t' 

:j(c  - l M*)  e 

li-l 


(2.4) 


where  tg  is  set  equal  to  0+,  the  time  at  which  the  delta  function  turns 
off.  Note  here  that  for  convenience  the  coupling  coefficient  and  the 
natural  mode  have  been  combined  into  one  term  A^(r).  The  s-plane  version 
of  Expression  (2.3)  is  written 


00  A.  (r) 

I(s,?)  = l rb~r  + Bi(?) 

i-1  S i 1 


(2.5) 


where  B^(r)  is  the  inverse  Laplace  transform  of  B^r)  6(t),  which  is  a 
constant  in  s and  a function  of  position  r. 

In  Equations  (2.1)  - (2.5)  the  series  contains  an  infinite  number 
of  terms.  In  general  only  the  first  few  terms  of  this  series  are  needed 
to  adequately  represent  the  late  time  response  of  the  system  [1],  [51 . 
The  early  time  response  on  the  other  hand  requires  a larger  number  of 
poles  for  reasonable  convergence  to  the  true  response.  This  is  in- 
tuitively reasonable  if  one  realizes  that,  in  general,  as  the  frequency 
or  the  imaginary  part  of  the  pole  increases  so  does  the  damping  constant 
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or  real  part.  Thus,  for  early  times  both  high  and  low  frequency  components 
are  present  and  as  time  progresses  tne  higher  frequency  components  damp 
out  and  disappear  until  only  the  first  few  dominant  resonant  modes  remain. 
Moreover,  all  transient  response  data  that  can  be  generated  either  ex- 
perimentally or  numerically  are  necessarily  band  limited  and  thus  contain 
only  a finite  number  of  poles.  For  these  reasons  the  series  is  truncated 
after  N terms,  where  N is  the  number  of  resonant  frequencies  contained 
in  the  transient  response  being  studied. 

2.2  Response  Due  to  an  Arbitrary  Excitation 

Since  a true  impulse  excitation  function  cannot  be  realized,  it  is 
of  interest  to  determine  the  form  of  the  transient  response  due  to  an 
arbitrary  exciting  waveform.  A general  response  function  R(s,r)  is 
given  in  the  Laplace  transform  domain  as 

R(s,r)  = F(s, r)  H(s,r)  (2.6) 


where  H(s,r)  and  F(s,r)  are  the  Laplace  transforms  of  the  system's  im- 
pulse response  and  the  arbitrary  driving  function,  respectively.  The 
system's  impulse  response  H(s,r)  was  written  in  (2.1)  and,  thus,  R(s,r) 
can  be  expressed  as 


R(s,r)  = F (s 


( N A (?) 


(2.7) 


where  the  possible  entire  function  has  been  neglected  and  the  coupling 
coefficient  and  natural  mode  have  been  combined  into  one  term,  A^(r). 

For  convenience  A^(r)  can  be  regarded  as  the  residue  for  the  ittl  complex 
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pole.  Note  that  the  residue  is  a function  of  position  on  the  object. 

If  the  inverse  Laplace  transform  of  (2.7)  is  performed,  then  the  tran- 
sient response  function  R(t,r)  is  written  in  general  form  as 
N _ s .t 

R(t,r)  = l B (r)  e 1 + g(t,r)  (2.8) 

i*l  1 

where 


B.(r)  = F(Si>?)  At(r) 


and  the  term  g(t,r)  is  dependent  on  the  form  of  the  driving  function. 

Equation  (2.8)  shows  that  the  transient  response  due  to  an  arbitrary 
exciting  waveform  can,  in  general,  be  written  as  a sum  of  complex  ex- 
ponentials plus  some  added  term  g(t,r).  The  desire  here  is  to  be  able 
to  express  the  transient  response  as  a sum  of  exponentials  only.  Thus, 
the  character  of  the  terra  g(t,r)  needs  to  be  studied  to  determine  if  it 
can  be  either  removed  or  expressed  also  as  a sum  of  exponentials. 

If  the  exciting  waveform  itself  has  pole  singularities,  as  in  the 
case  of  a step  function  or  a sinusoidal  function,  then  g(t,r)  can  be 
expressed  as  a sum  of  exponentials.  As  an  example,  consider  a driving 
function  which  is  a simple  step  function  u(t),  then,  the  Laplace  trans- 
form of  the  driving  function  is 
F(s)  = 1/s 


If  the  impulse  response  of  the  system  is 
N A . 


H(s)  * l 


i-1  S “ Si 


then  the  response  R(s)  is  simply 

1 N A< 

R(s)  - F (s)  H(s)  = - [ — 

Si«l  S " Si 


(2.9) 


8 


— - ■■  ■ _ 


When  the  inverse  Laplace  transform  is  performed  on  R(s),  one  obtains 


a-1°°  , N A, 

R(t)  - ■—  / - l — eSt  ds 

2nj  aijco  siti  3 - si 


u(t) 


N 

l \ 

i=l 


e 1 


Si  SiJ 


u(t) 


N s .t 

I Bie 

i=0 


where 


Bo  • -Jj  Ai/sl 


and 


(2.10a) 


(2.10b) 


(2.10c) 


Thus,  since  the  step  driving  function  has  a singularity  at  the  origin 
then  the  response  function  is  written  as  a simple  sum  of  exponentials. 

If  the  exciting  waveform  is  of  finite  duration,  that  is,  it  is 
turned  off  after  some  time  tg,  then  it  can  be  shown  that  the  term  g(t,r) 
is  identically  equal  to  zero  for  time  t greater  than  tg.  Consider,  for 
example,  the  case  where  the  driving  function  is  a simple  pulse 


Again  let  the  impulse  response  be 


H(s)  = l 


i-1  S ' Si 


then,  the  response  function  in  the  Laplace 


1 stn  N 

R(s)  = F(s)  H(s)  - ±(1  - e U)  £ 


i=l 


transform  domain  is 

A. 

1 


(2.11) 


The  corresponding  time  domain  expression  for  the  response  function  is 


R(t)  « 

o+j® 

/ R(s)  est 

ds 

(2.12a) 

Q 

1 

L-u 

8 

A 

s.t 

R(t)  -< 

y -i 

u s 
i=l  Si 

(e  1 - 

o 

XI 

V 
XI 

V 

o 

H 

(2.12b) 

N 

S t 

l Bi 

li-1 

X 

e 

rr 

V 

rr 

O 

where 


~Si\ 

e ) 


Expression  (2.12b)  shows  that  after  time  tg  the  response  function  is 
simply  a sum  of  the  exponentials  which  characterize  the  body  while,  be- 
fore tg,  the  sum  of  exponentials  contains  one  exponential  term  which  is 
dependent  on  the  driving  function. 

When  the  driving  function  is  not  finite  in  time  and  has  no  pole 
singularities,  as  in  the  case  of  the  Gaussian  pulse,  then  the  g(t,r) 
term  never  disappears  and  cannot  be  written  as  a sum  of  exponentials. 
However,  this  difficulty  may  be  circumvented  by  simply  deconvolving  the 


response  function  R(t,r)  in  a standard  manner  to  obtain  the  system's 
impulse  response  H(t,r).  Equation  (2.6)  gives  for  instance 

H(S)  = F(ij  * (2.13) 

Thus,  removing  the  driving  function  from  the  response  function  results 
in  the  impulse  response  which  is  assumed  to  be  a sum  of  exponentials 
only.  The  inevitable  presence  of  experimental  and  computational  noise 
limits  the  upper  frequency  for  which  the  deconvolved  spectrum  H(s)  is 
accurate  and,  in  practice,  the  computed  spectrum  must  be  truncated  be- 
yond this  frequency.  This  simply  sets  a limit  on  N.  Care  also  must 
be  taken  to  exclude  the  time  t = 0 from  the  deconvolution  because  of 
the  presence  of  the  delta  function  discussed  previously. 

Numerical  examples  of  the  above  cases  are  given  in  Chapter  3 once 

V 

the  method  by  which  the  poles  will  be  extracted  has  been  developed. 
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3.  THE  NUMERICAL  METHOD  FOR  EXTRACTING  POLES 


* 


I 


In  the  previous  chapter  it  was  shown  that  the  impulse  response  and 
for  certain  cases  the  general  transient  response  of  a system  can  be  ex- 
pressed as  a sum  of  a finite  number  of  exponentials.  This  chapter  will 
develop  the  numerical  method  by  which  the  values  of  the  poles  and  their 
corresponding  residues  can  be  obtained  if  the  impulse  or  transient  re- 
sponse of  the  system  is  given.  The  numerical  method  used  is  based  on 
Prony's  algorithm  [12],  [13],  [14]  which  does  not  seem  to  be  widely  known. 
Prony's  algorithm  has  been  used  in  the  fields  of  automatic  control  [16] 
and  biological  signal  processing  [17],  but  only  to  represent  or  synthesize 
a signal  in  terms  of  a set  of  exponentials  which  do  not  necessarily  have 
any  physical  relationship  to  the  system  which  produced  the  signal.  The 
desire  here,  however,  is  to  extract  from  a system's  transient  response  a 
set  of  complex  exponentials  which  are  in  fact  the  characteristic  reso- 
nances of  the  system  being  studied. 

In  this  chapter  Prony's  method  is  developed  for  a system  with  simple 
poles  only.  Prony's  method  is  then  extended  to  systems  containing  multi- 
ple poles.  The  derivation  of  Prony's  method  is  presented  in  detail  in 
several  different  ways  so  that  the  reader  will  get  a good  feeling  for 
the  method  and  the  problems  associated  with  it.  An  abbreviated  treatment 
is  presented  in  [18].  Several  numerical  examples  are  presented  in  order 
to  establish  some  guidelines  for  the  use  of  the  method  and  to  illustrate 
some  of  the  problems  associated  with  the  method. 


! 
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* 


a 


12 


ET 


A^  + A2  + . . . 

+ ^ 

A^Z^  + ^2^2 

+ Vn 

A1Z12  + A2Z2  + ‘ “ * 

+ ANZN 

Vl  - A1Z1  " + A2Z2  Vn 


M-l 


(3.3a) 


where 


Zi  = 6 


s^t 


(3.3b) 


In  this  set  of  equations  it  is  necessary  to  solve  for  both  the  N values 
of  Z^  and  the  N values  of  the  A^.  This  solution  requires  that  the  value 
of  M be  at  least  equal  to  2N  ; however,  the  solution  to  this  set  of 
equations  is  nontrivial  since  they  are  nonlinear  in  the  Z^'s. 

Let  Z^,  . . . , ZN  be  the  roots  of  the  algebraic  equation 


+ ci^Z^  + + . . . + ct^Z^  * 0 


(3.4a) 


so  that  the  left-hand  side  of  (3.4)  is  equal  to  the  product 


(Z  - zL)  (Z  - z2)  . . . (Z  - zN) 


that  is , 


N n 

l “„z"  -i2i  <z  - V 


m»0 


(3.4b) 


The  coefficients  may  be  determined  as  follows.  Multiply  the  first 
equation  in  (3.3)  by  a^,  the  second  by  . . .,  the  N + 1 by  and 

add  the  results.  Since  each  of  the  satisfies  Equation  (3.4a),  then 


the  result  is  of  the  form 


aoRo  + alRl  + • • • + Vn  3 0 • 


A set  of  M - N - 1 additional  equations  is  obtained  in  the  same  manner 
by  starting  successively  with  the  second,  third,  . . .,  (M  - N)Ctl 
equations.  Thus,  it  is  possible  to  obtain  the  M - N linear  difference 
equations 


a0R0  + alRl  + • • • “A  = ° 


(XqR^  + =^R2  + • • 


“n^i+i 


a0RM-N-l  + + 


aNRM-] 


which  can  be  more  simply  written  as 


I aR.,=0;p+k=n=  0,1,  . . .,  M - 1 
P P+k 


Thus,  the  sampled  values  R^  satisfy  an  N order  linear  difference  equation. 
This  difference  equation  will  be  referred  to  as  Prony's  difference  equation. 


In  Equations  (3.5)  the  R^  are  known,  hence,  this  set  can  be  solved 


to  obtain  the  N + 1 coefficients,  a.  If  M = 2N  then  the  system  of  equa- 
tions may  be  solved  directly  by  matrix  inversion,  and  if  M > 2N  the  system 
may  be  solved  approximately  by  the  method  of  least  squares.  Equations 


(3.5)  are  most  conveniently  solved  by  defining  at^  equal  to  one;  then, 


the  ap's  may  be  obtained  by  solving  the  equations 


N-l 


l Vp+k  = _EW;  P + k = n = °*1-  • • •’  M - 1 

p=0  r 


(3.6) 


If  M = 2N  this  set  of  equations  is  written  in  matrix  form  as 


AB  = C 


where 


R0  R1  R2 


R1  R2  R3 


R2  R3  R4 


R-  R.  R_ 

3 4 3 


^-1  ^ 


^-1 


. . R 


2N-2 


16 


(3.7a) 


(3.7b) 


Matrix  A is  a square  symmetric  circulant  matrix  and  thus  is  readily  in- 
vertible. If  M is  greater  than  2N  , then  the  set  of  Equations  (3.6) 
can  be  solved  using  a least-squares  approach.  This  is  most  conveniently 
done  by  performing  a pseudo-inverse  to  write  Equations  (3.6)  in  matrix 
form  as 

AT  AB  » AT  C (3.3a) 

where  A is  now  a rectangular  matrix  of  the  form 
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*n 

Vt-i 


IV 


(3.8c) 


The  square  matrix  formed  by  multiplying  matrix  A by  its  transpose  A is 
simply  a Gramian  matrix  and  is  also  symmetric  and  circulant. 

Once  the  a's  have  been  determined  by  either  of  the  above  approaches, 


the  next  step  is  to  solve  for  the  N values  of  Z. . These  Z.  are  ob- 

i i 


tained  by  finding  the  roots  of  the  polynomial  (3.4a).  The  N roots  are 
complex  numbers  and  appear  in  complex  conjugate  pairs.  The  roots  of 
(3.4a)  may  be  found  by  using  any  of  the  several  polynomial  root-finding 
methods.  The  most  accurate  and  quickest  of  these  methods  is  a routine 
based  on  Muller's  method  [19],  [20]. 

It  is  now  a trivial  matter  to  obtain  the  poles  s^.  Since  the"Tbots 
of  (3.4a)  were  defined  by  (3.3b)  as 


.! 
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then  the  poles  are  simply 


(3.9) 


The  final  step  in  Prony's  method  is  to  determine  the  values  of  the 
residues  A^.  To  do  this  one  simply  solves  the  matrix  equation  embodied 
in  Equations  (3.3a).  In  matrix  form  this  set  of  equations  is  written  as 


DE  - F 


inhere 


1 1 


D - Zt  Z2  • • • Z„ 


2 2 2 
Z1  Z2  ‘ ’ ZN 


_N-1  _N-1 

^1 


E - A„ 


(3.10c) 


Vi  J 


(3. lOd) 


Note  that  the  matrix  D is  in  the  form  of  a Vandermonde  matrix  which  has 
an  inverse  which  can  be  computed  in  closed  form.  If,  however,  one  wishes 
to  solve  Equation  (3.3a)  by  using  a least-squares  approach  then  the  re- 
sultant matrix  equation  has  a matrix  which  is  only  symmetric  and  cannot 
be  inverted  as  simply  as  the  Vandermonde  matrix  D. 


3.1.2  Relationship  to  the Z -transform 

The  method  presented  in  the  above  subsection  for  finding  the  simple 
poles  and  residues  from  a set  of  discrete  transient  data  is  Prony's 
method  as  it  is  usually  derived  and  used.  However,  Weiss  and  McDonough 
[21]  have  demonstrated  that  Prony's  method  may  be  regarded  as  a Pade 
approximation  in  the  Z-transform  domain. 

Consider  the  Z-transform  of  R(t)  to  be 

R(Z)  = Rq  + RjZ"1  + R2Z~2  + . . . + R2N-lZ~(2N~i:>  + (3.11) 


In  the  Pad£  approximation  one  seeks  to  equate  the  first  2N  terms  of  (3.11) 
with  a function  of  the  form 


F (Z) 


V + Viz 


N-l 


. + a2Z 


,N  _N-1 

Z + VlZ  + 


. + a^Z  + Oq 


(3.12) 


20 


Note  that  the  denominator  in  (3.12)  is  the  same  as  the  polynomial  in 
(3.4a).  Carrying  out  the  Padd  procedure  by  equating  (3.11)  and  (3.12) 
and  multiplying  through  by  the  denominator  of  (3.12)  yield 


a^Z  + a^^Z1  + . . . + a^Z  “ (Z‘  + <*jj_jZ  + . . . + a^Z  + ciq 


(Rq  + R^"1  + R2Z-2  + . . . + r2n_1z~(2N_1)  + • * ■*> 


Equating  the  coefficients  of  like  powers  of  Z in  (3.13)  yields  the  fol 
lowing  set  of  equations. 


^-1  - R0  aN-l  + R1 

*N-2  = R0  “N-2  + R1  “n  + R2  > 

al  = R0  “1  + R1  a2  + • ' * + %-l  > 

Ro  ao  + R1°1  + • • • ^-l  Vi  + * 0 

R1  a0  + R2  al  + * ' * S V-l  + Vn  = 0 

> 

^-1  a0  + S °1  + ■ * * R2N-2  “n-I  + R2N-1  ’ ° 


) 

(3.13) 


(3.14a) 


(3.14b) 
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Equations  (3.14b)  are  the  same  as  Equations  (3.5a)  of  the  previous 


L 


section. 

If  the  Z-transform  of  Equation  (3.1)  is  taken,  the  result  is 


N 


R(Z)  = l A 
i-1 


Z 


z - z. 


(3.15) 


where  Z^  = e 


s^nAt 


Hence,  the  and  the  Z^  of  the  Pad4  approximation 
in  the  Z-transform  domain  are  the  same  as  the  and  the  Z of  the  pre- 
vious section.  This  method  allows  one  to  solve  for  the  by  a somewhat 
simpler  method.  Rather  than  solving  the  matrix  Equation  (3.3a),  one 


can  simply  solve  for  the  a^  in  (3.14a)  and  perform  the  partial  fraction 


R (Z  ) 

expansion  of  ^ ' to  yield 


1 N 

z -l  T 


(3.16) 


3.1.3  Relationship  to  Corrington's  difference  equation 
In  1965  Corrington  derived  a difference  equation  from  which  it  is 
possible  to  extrapolate  the  time  response  of  a linear  lumped-constant 
time-invariant  network  to  late  time  by  knowing  only  discrete  values  of 
the  early-time  response  [22].  His  difference  equation  can  be  shown  to 
be  Prony's  difference  equation  in  a somewhat  camouflaged  form.  His 
derivation  is  interesting  and  is  presented  here  in  a similar  form. 

This  derivation  is  useful  when  the  development  of  the  derivation  for 
multiple  poles  is  presented  in  the  following  section. 

Equation  (3.2)  may  be  rewritten  in  an  alternate  form  as 


N 


R[ (n  - r)  At] 


s.nAt  -s  rAt 

7 A|  e e ; r - 0,1, 


N 


(3.17) 


i-1 
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a-  _ 


The  term  A^e  * can  be  eliminated  from  the  above  N + 1 linear  equations 

by  subtracting  each  equation  in  turn  from  the  preceding  one  after  multi- 
set 

plication  by  e . This  gives  the  set  of  N linear  equations 


R[  (n  - r)  At)  - e R[  (n  - r - 1)  At] 


s,At  N-l  -rAts.^.  i / -s,  At  -s,  Af 
J 1 l e 1+1  ± e k -e  k+1 

i-1  k=1  V 


s nAt 

• A.+1  e 1 x ; r - 0,1 N-l 

s2n^t 

from  which  the  term  A2  e can  be  eliminated  by  the  same  procedure. 

The  result  is 


/ s At  s,,At\ 

R((n  - r)  At]  - R[ (n  - r - 1)  At]  fe  + e 1 + R[ (n  - r - 2)  At 


s^At  s2At \ s2At  N-2  -rAtsi+2  2 
i e * e T e ,11, 

' i-1  k-1 


I “skAt  -Sk+2AM  A Si+2nAt 

I • A1+2  e 


r*  0,1,2,  . . .,N-2 


s^nAt 


If  this  process  is  continued  until  all  N of  the  A^  e are 


eliminated,  the  result  is 


l Vp  R([n  " Pl  At) 

p-0  y 


where  the  N + 1 values  of  a are  the  exact  same  a appearing  in  the 

P P 

algebraic  Equation  (3.4b).  Since  the  Z.  are  defined  as 


s .At 

then  the  are  sums  and  products  of  the  e . For  example  if  N = 3, 
then  (3.4b)  is  written 


s,  At 


l . Z»-  Z-.'1  ||z-.2  Hz-.2 

p=0  p 


s„At 


s„At 


and,  thus  the  a's  are 


°3  = 1 

s^At  S2At  s At 
= -e  - e - e 


s^At  S2At  s^At  s^At  S2At  s.At 
a^“e  e +e  e +e  e 


s^At  s2At  s^At 
«g  = -e  e e 


Hence,  Equation  (3.20)  is  an  alternative  way  of  writing  Prony's  difference 
Equation  (3.5b).  Corrington  writes  the  difference  equation  as 
N 

- ijj  R(nAt)  = 7 aN_  R[(n  - p)  At];  n >_  N (3.21) 

P=1  '~P 

which,  if  the  a's  are  known,  expresses  the  repsonse  at  some  time  n in 
terms  of  N previous  time  sampled  values  of  the  response. 


3.2  Prony's  Method  for  Multiple  Poles 

For  the  most  part  electromagnetic  antennas  and  scatterers  possess 
only  simple  poles.  However,  Tesche  [23]  has  shown  that  a dipole  can  be 
resistively  loaded  in  such  a way  as  to  make  it  critically  damped,  that 
is,  to  have  a double  pole  on  the  negative  real  axis.  Multiple  poles 
also  result  in  the  transient  response  of  a system  if  the  system  is 
driven  by  a signal  which  itself  has  a multiple  pole.  The  most  common 
example  of  this  would  be  a system  excited  by  a ramp  waveform.  The 
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ramp  waveform  has  a double  pole  located  at  the  origin.  It  is  shown  in 
this  section  that  Prony's  method  can  be  used  to  extract  these  multiple 
poles  without  knowing  £ priori  that  they  are  present. 

The  general  form  of  the  transient  response  for  a system  containing 


both  multiple  and  simple  poles  can  be  written  as 

N ( Mi  s.t 

R(t)  = l <1  + l tJ'  B P(M  ) } A e 

i=l  V.  j=2  J J 1 


(3.22) 


where 


P(Mi)  = 0, 

if  Mi  < 2 

P(Mi)  = 1, 

if  M £ 2 

and  where  is  the  multiplicity  of  the  it^1  pole.  For  example,  if  the 
iCh  pole  is  a double  pole,  then  = 2 and  P(frL)  = 1.  In  discrete  data 

form.  Equation  (3.22)  can  be  written 

M, 


N i s 

R[(n  - r)  it]  ■ W 1 + [ [ (n  - r)  At ]J_X  B P(M.)>  A.  e 1 

i=l  1 ' j=2  J 


s (n-r)At 


r = 0,1,2,  . . . , L 

(3.23) 


where  L is  the  total  number  of  poles  if  each  pole  is  taken  times. 
If  as  for  the  simple  pole  case 


Zi-* 


s^  At 


then  an  algebraic  equation  can  be  written 

N M.  L 

n (Z  - Z.)  1 = Y a zm  - 0 . 

1 _ 1 m 


i*l 


m^O 


(3.24) 
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It  is  shown  that  the  L Equations  (3.23)  can  be  reduced  to  give  Prony's 
difference  equation 


£ aT  _ R[  (n  - m)  At]  » 0 


(3.25) 


This  difference  equation  is  derived  in  the  same  manner  as  in  Section  3.1.3 

except  for  one  small  difference.  As  before,  each  equation  in  turn  is 

s . At 

subtracted  from  the  previous  one  after  multiplying  it  by  e . This 

Si+lAt 

"tep  is  repeated  bh  times  and  then  e is  used.  This  can  be  best 

demonstrated  in  the  following  example. 

Let  the  response  be 

2 f Mi  , . ^ s.(n-r)  At 

R[(n  - r)  At]  = £ ( 1 + £ [(n  - r)  At]J_i  B P(M  )}  A e ; 

i-1  j>2  J1  M 1 

r - 0,1, 2, 3 


and  if  At  ■ 1,  M * 2,  = 1,  then, 


s (n-r)  s (n-r) 

R[n  - r]  = (A  + Bt)  e + C e 


r = 0,1, 2, 3 


(3.26) 


Note  here  that  the  constant  B is  simply  the  product  B_,^A^.  If  one  uses 
the  method  of  Section  3.1.3,  the  four  Equations  (3.26)  are  reduced  to 
the  three 


s1  s n -rs  s.n  -rs? 

R(n  - r)  - e R(n  - r - 1)  ° B e e +Ce  e 

s -s 

• (1  - e 1 e “)  , 
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S1  Sin  ~rS1  _S1  S->n  ~rS9 

R(n  - r - 1)  - e 1 R(n  -r-2)  = Bei  e e 1 + C e 2 e 2 


“S2  "2s2  S1 
( e 2 - e 2 e X)  , 


S1  sin  “rsi  ~2si 

R(n  - r - 2)  - e A R(n  -r-3)  = BeX  e e L 


s n -rs,  -2s  -3s,  s. 

+ C e ^ e • (e  - e d1), 


(3.27) 


Since  the  pole  s^  was  a double  pole,  the  step  is  repeated  by  multiplying 
S1 

by  e again.  The  two  resulting  equations  are 

S1  S1  2sl 

R(n  - r)  - e R(n  - r - 1)  - e R(n  - r - 1)  + e R(n  - r - 2) 


s_n  -rs, 
- C e e 


, S2  S1  "2s2  2s1. 
(l-2e  e +e  e ) 


(3.28a) 


S1  S1  2sl 

R(n  - r - 1)  - e R(n  - r - 2)  - e 1 R(n  - r - 2)  + e 1 R(n  - r - 3) 


-2s,,  s - -3s,,  2s 

C e “ e ‘ (e  ‘ - 2 e ^e+e  e).  (3.28b) 


s,n  _rs2  — s , 


Finally,  the  difference  equation  is  obtained  by  multiplying  (3.28b)  by 
S2 

e and  subtracting  from  (3.28a).  The  difference  equation  is  thus 


S-  s , s. 


2s, 


R(n  - r)  - (e  1 + e 1 + e 2)  R(n  - r - 1)  + (e  1 + 2 e 1 2)  R(n  - r - 2) 

(3.29) 


2s  s 

- (e  2 e 2)  R(n  - r - 3)  - 0 


From  (3.29)  the  a's  of  (3.25)  are 


°3  ” 1 


S1  S1  s2 
e - e - e 


2sl  S1S2 

a1-e  1 + 2 e X 

2S1  S2 
°0  " -®  e 
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This  example  shows  that  although  there  were  only  two  distinct  poles  the 
difference  equation  for  the  response  was  of  order  three  producing  four 
coefficients. 

In  the  above  example,  it  was  shown  that  Prony's  difference  equation 
is  also  valid  for  transient  responses  containing  multiple  poles.  The 
order  of  the  difference  equation  does  increase  from  N to  L thus  requiring 
L + 1 values  of  a to  be  solved  for.  Once  the  L + 1 values  of  a have 
been  determined,  then  the  next  step  is  to  find  the  L roots  of  the  poly- 
nomial (3.24).  If  Muller's  method  is  used  for  finding  the  L roots  of 
the  polynomial,  it  does  indeed  return  an  individual  root  as  many  times 
as  its  multiplicity  requires.  Therefore,  after  finding  the  roots  of 
the  polynomial,  the  roots  are  scanned  to  see  if  any  appear  more  than 
once,  indicating  the  presence  of  a multiple  pole.  Hence,  it  is  possible 
to  proceed  with  Prony's  method  to  the  point  of  finding  the  poles  with- 
out knowing  if  there  are  multiple  poles  present  or  not.  After  the  poles 
and  their  multiplicity  have  been  determined,  then  the  problem  is  to 
calculate  the  values  of  the  residues. 

The  calculation  of  the  residues  is  done  by  solving  the  equation 

N f Mi  , . ) s.nAt 

R(nAt)  = l < 1 + £ (nAt)J_1  B P (M  ) ) A e 1 n - 0,1,  . . . (3.30) 
i=l  ( j=2  J j 

which  differs  from  Equation  (3.2)  because  of  the  presence  of  the  terms 
(nAt)^  \ This  fact  causes  the  matrix  Equation  (3.10)  to  be  changed  by 
multiplying  certain  columns  by  multiples  of  (nAt).  For  example,  if 
from  the  previous  example. 
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s.  At 

Z1  ' 6 

s.  At 

Z2  - e 

s_At 

Z3  " 6 

then  the  resulting  matrix  equation  to  be  solved  in  order  to  find  the 
coefficients  A,  B,  and  C of  (3.26)  is 


f — 

10  1 

A 

R(0At) 

Z1  lAtZ2  Z3 

B 

= 

R(lAt) 

2 2 2 
Zx  2AtZ2  Z3 

C 

R(2At) 

w 

■ — — • 

^ _ 

where  column  2 has  been  changed  because  of  the  presence  of  the  term  Bt 
in  (3.26).  The  next  example  shows  that  this  necessity  to  alter  certain 
columns  of  the  matrix  leads  to  instabilities  in  the  solution  process. 

As  a more  complicated  example  of  a system  with  a double  pole  con- 
sider the  transient  response 

R(t)  = 3.0  + 7.0  t + 6.0  e cos(4t  - w/6)  + 4.0  e cos(6t  + n/3) 

(3.31) 

which  is  plotted  in  Figure  3.1.  Note  that  the  second  term  in  (3.31)  is 
a ramp  response  which  is  very  dominant  in  the  plotted  transient  response. 
This  ramp  term  gives  rise  to  a double  pole  at  the  origin.  When  Prony's 
method  was  applied  to  these  data,  the  extracted  poles  and  residues  in 
Table  3.1  were  obtained.  Note  that  the  extracted  poles  are  all  within 
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SAMPLE  VALUES 
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TIME  SAMPLES  (At  = 0.1  sec.) 

Figure  3.1.  Response  containing  a double  pole  at  the  origin. 
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EXTRACTED  POLES  AND  RESIDUES  FOR  THE  DOUBLE  POLE  CASE 


0.08  percent  of  the  original  poles  and  that  the  double  pole  was  indeed 
found.  However,  the  residues  that  were  recovered  do  not  compare  with  the 
true  residues  at  all.  The  residues  do  not  even  fall  in  complex  conjugate 
pairs  as  is  required.  This  is  undoubtedly  due  to  the  fact  that  the  first 
matrix  was  altered,  as  discussed  in  the  previous  example,  in  such  a way 
as  to  make  it  ill-conditioned.  The  important  point  in  this  example  is, 
however,  that  the  double  pole  was  extracted  from  the  transient  response. 


3.3  Numerical  Examples 

In  using  Prony's  method,  as  in  the  use  of  any  numerical  routine, 
several  guidelines  should  be  followed  in  order  for  the  method  to  work 
accurately  and  quickly.  For  the  most  part,  these  guidelines  were  ob- 
tained after  running  several  examples  on  the  computer  and  studying  the 
results.  Therefore,  this  section  presents  several  numerical  examples 
which  point  out  the  problems  and  guidelines  which  one  must  follow  for 
the  successful  use  of  Prony's  method. 

In  this  section,  for  all  but  one  example,  the  data  analyzed  are 
from  the  transient  response  of  the  current  on  a 1.0  m long  dipole  with 
a half-length-to-radius  ratio  of  100.  The  data  were  obtained  using  a 
numerical  time-domain  computer  code  [24].  The  antenna  was  modeled 
using  sixty  equal-length  segments  and  the  exciting  field  was  a Gaussian 

pulse  which  was  applied  across  the  center  two  segments  of  the  antenna 

2 2 

model.  The  Gaussian-pulse-time  variation  was  exp (-a  (t  - t ) ) with 

r max 

9 -1 

a,  the  Gaussian  spread  parameter,  equal  to  5 x 10  s and  t^^  equal 
to  5.556  x 10  ^ s.  The  induced  current  at  the  center  of  each  of  the 
sixty  segments  was  calculated  for  500  time  steps,  where  the  time  step 
size  was  At  = 5.556  x 10  ^ s.  Note  that 


L 

60  x C 
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where  L is  the  length  of  the  structure,  1 in,  and  C is  the  speed  of  light. 
In  order  to  calculate  the  complex  resonances  of  the  structure,  the  cur- 
rent on  one  of  the  center  or  source  segments  is  used.  This  current  is 
plotted  in  Figure  3.2. 

The  study  of  Figure  3.2  reveals  several  facts  about  the  dipole  re- 
sponse. The  most  apparent  is  the  damped  oscillatory  behavior  of  the  cur- 
rent. The  initial  spike,  which  is  a maximum  at  t = 10At,  is  the 

max 

Gaussian  driving  function  being  applied.  This  spike  is  followed  by  an 
immediate  negative  spike  which  decays  toward  zero  until  time  equals  60At. 
The  response  from  time  0 to  60At  is  very  similar  to  the  response  of  an 
infinite  long  dipole.  This  response  is  expected  since  the  current  which 
propagates  down  the  antenna  and  gets  reflected  back  from  the  end  is  not 
seen  at  the  center  of  the  antenna  until  60At.  Thus,  the  current  at  the 
center  of  the  antenna  is  not  affected  by  the  length  of  the  antenna  until 
the  current  has  propagated  to  the  end  and  back. 

Example  1 

The  first  example  of  the  application  of  Prony's  method  is  shown  in 
Figure  3.3.  Of  the  500  current  samples  available,  only  eighty  sampled 
values  were  actually  used.  These  eighty  samples  were  taken  from  the 
first  160  current  samples  at  every  second  time  step.  Figure  3.3a  shows 
the  range  of  the  current  which  was  used  to  fill  the  matrix  of  Prony's 
difference  equation.  Prony's  method  was  solved  using  the  standard  in- 
version technique.  Thus,  since  eighty  sampled  values  were  used,  forty 
poles  and  residues  were  obtained.  Figure  3.3b  shows  the  left  half  of 
the  s-plane  in  which  the  forty  extracted  poles  have  been  plotted. 

Figure  3.4  shows  a comparison  of  the  extracted  poles  with  the  true  even 
poles  of  the  dipole  as  calculated  by  Tesche  [5].  As  can  be  seen  only 
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Figure  3.2.  Source  current  on  1.0  m dipole. 
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Figure  3.3.  (a)  Data  window  used  in  Example  1. 

(b)  Locations  of  extracted  poles  in  Example  1 
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eleven  pole  pairs  correspond  to  the  true  poles  of  the  system  leaving 
eighteen  poles  which  have  no  relationship  to  the  system.  The  fact  that 
only  the  first  ten  even  poles  were  generated  is  not  surprising  if  one 
considers  the  original  model.  The  original  model  was  a thin-wire 
approximation  and  was  driven  with  even  symmetry  by  a Gaussian  pulse. 
Because  of  the  thin-wire  model  and  the  width  of  the  Gaussian  pulse,  the 
expected  spectral  response  has  an  upper  frequency  limit  of  about  10  L/X 
[25],  where  L equals  the  antenna  length.  The  resonances  for  a dipole 
occur  at  approximately  X/L  ■ 1/2,  3/2,  5/2,  ...»  thus,  with  the  upper 
frequency  limit  of  10  L/X  not  many  more  than  the  first  ten  even  reso- 
nances can  be  expected.  The  extra  eighteen  poles  appear  for  two  reasons. 
In  Prouy's  method,  if  the  least  squares  approach  is  not  used,  when  2N 
sampled  values  are  used  the  method  returns  N poles.  Thus,  Prony's 
method  is  forced  to  return  more  poles  than  are  present  in  the  system. 
Also,  since  the  Gaussian-pulse  driving  function  is  only  a model  of  a 
delta  function, then  the  response  function  is  not  a true  impulse  response 
but  is  a response  of  the  form  of  (2.8).  Since  the  g(t)  term  of  (2.8) 
was  not  removed,  the  extra  poles  present  are  needed  to  represent  it. 
Figure  3.5  demonstrates  how  well  the  transient  response  was  modeled  with 
these  forty  poles  and  their  associated  residues.  Note  that  1000  time 
steps  were  generated  when  only  the  first  160  of  the  original  time  steps 
were  used  to  generate  the  poles.  This  shows  how  Prony's  method  can  be 
used  to  extrapolate  late  time  data  from  a small  set  of  early  time  data. 
Figure  3.6  is  a three-dimensional  plot  of  the  second  quadrant  of  the 
s-plane  showing  both  the  position  and  the  amplitude  of  the  poles.  The 
length  of  the  vertical  lines  represents  the  amplitude  of  the  poles  on 
a logarithmic  scale.  The  amplitude  of  the  true  poles  tapers  off  as  the 
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Figure  3.5.  Reconstructed  transient  response  using  extracted  poles 
of  Example  1. 
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Figure  3.6.  Position  anti  amplitude  of  poles  in  the  second  quadrant 
of  the  s-plane  for  Example  1. 
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frequency  increases  because  the  amplitude  of  the  Gaussian  driving  function 
in  the  spectral  domain  also  tapers  off  as  frequency  increases.  It  should 
be  noted  that  the  amplitudes  of  some  of  the  extra  poles  are  equal  to  or 
greater  in  magnitude  than  those  of  the  true  poles.  These  extra  poles  with 
the  large  amplitude  are  the  poles  which  are  predominantly  representing 
the  driving  function. 

Example  2 

As  a second  example  only  the  first  sixty  samples  of  the  response 
are  used.  This  sampling  interval  is  only  large  enough  to  contain  the 
response  resembling  the  infinite  wire,  Figure  3.7a.  The  resulting  poles 
are  plotted  in  Figure  3.7b.  Note  that  none  of  these  poles  are  the  true 
poles  of  the  system  but  that  the  sampled  portion  of  the  response  is 
accurately  reproduced  using  these  poles,  Figure  3.8.  Also,  note  that 
only  twenty-five  poles  are  plotted  in  Figure  3.7b  when  thirty  were  ex- 
pected. This  is  due  to  the  fact  that  the  program  was  written  so  as  to 
eliminate  any  right-half-plane  poles  that  appear.  In  this  case,  five 
poles  with  positive  real  parts  were  removed  before  the  residues  were 
calculated.  This  example  points  out  the  necessity  to  include  some  por- 
tion of  the  oscillatory  part  of  the  response.  It  was  found  that  at 
least  one  cycle  of  the  lowest  frequency  present  must  be  used.  Thus, 
the  sample  interval  of  example  1 is  the  minimum  which  can  be  used  to 
obtain  the  true  poles  of  the  system. 

Example  3 

This  example  also  takes  sixty  samples  but  the  time  step  size  is  now 
3At  thus  giving  a sampling  interval  of  the  required  size.  Figure  3.9a, 
and  requiring  the  return  of  only  thirty  poles  as  opposed  to  forty  in  the 
first  example.  The  twenty-eight  poles  that  were  found  are  plotted  in 
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Figure  3.7.  (a)  Data  window  used  in  Example  2. 

(b)  Locations  of  extracted  poles  in  Example  2. 
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Figure  3.9b.  Note  that  the  true  poles  appear  while  the  extraneous  poles 
have  been  reduced  in  number.  The  extra  poles  still  must  appear  because 
of  the  presence  of  the  driving  function  in  the  transient  response. 

Figure  3.10  points  out  another  problem  inherent  in  the  solution  of  the 
matrix  equations  without  the  use  of  a least  squares  technique.  Since 
twenty-eight  poles  were  found,  it  was  necessary  to  solve  for  twenty- 

I 

eight  residues  using  the  matrix  Equation  (3.10).  Only- theT' first  twenty- 
eight  samples  of  the  current -wertfneeded  to  fill  the  vector  F of  (3.10a). 
Thus,  the  resulting  response  is  expected  to  be  accurate  at  only  the  first 
twenty-eight  time  steps.  This  is  represented  in  the  plot  of  Figure  3.10. 
The  obvious  correction  to  this  problem  is  to  use  more  samples  to  fill 
the  vector  F by  using  the  least  squares  method.  An  example  of  this 
approach  will  be  presented  later  in  this  section. 

Example  4 

The  Nyquist  criterion  [26]  states  that  the  sampling  rate  of  a tran- 
sient response  must  be  at  least  twice  as  fast  as  the  highest  frequency 
present  in  the  data  in  order  to  be  able  to  resolve  that  frequency.  That 
is,  the  sampling  rate  At^  is 


AtN  - 2 L 


where  |(  is  the  highest  frequency  desired.  In  the  previous  example  the 
time  step  size  was  3At  3 1.6668  x 10  ^ s.  Hence,  the  highest  frequency 
which  could  be  expected  was  L,  * 2.9998  x 10^  Hz  or  — ■ 19.9984  which 
is  indeed  just  higher  than  the  highest  frequency  pole  obtained.  Figures 
3.11  and  3.12  also  demonstrate  this  phenomenon.  In  Figure  3.11a  the 
time  step  size  is  4At  = 2.2224  x 10  ^ s and  the  resulting  frequency 
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Figure  3.10.  Original  response  minus  the  reproduced  response 
of  Example  3. 
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Figure  3.11.  (a)  Data  window  used  when  the  time  step  is  4At. 

(b)  Locations  of  the  resulting  poles. 
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limit  is  wL/cir  * 14.99.  The  highest  pole  frequency  obtained  is  14.80. 

Likewise  in  Figure  3.12a  the  time  step  is  5At  = 2.778  x 10  ^ s giving 
a frequency  limit  of  11.999  and  the  highest  frequency  pole  obtained  is 
11.4.  Thus,  in  order  to  obtain  the  highest  frequency  pole  that  is  con- 
tained in  the  data  one  must  sample  the  transient  response  at  at  least 
the  Nyquist  rate. 

B 

In  all  of  the  previous  examples  no  attempt  was  made  to  remove  the 
influence  of  the  driving  function,  the  Gaussian  pulse,  from  the  response 
function.  Thus,  poles  have  been  present  which  were  required  to  fit  the 
g(t)  term  of  Expression  (2.8).  Since  the  driving  function  is  a Gaussian 
pulse,  it  cannot  be  modeled  exactly  with  a finite  number  of  poles  but 
can  be  assumed  to  go  to  approximately  zero  for  some  value  of  time  tg. 

Hence,  if  the  response  function  is  sampled  after  time  tg,  the  g(t)  term 
should  be  zero.  Also,  since  the  analytical  form  of  the  spectrum  of  the 
Gaussian  pulse  is  known,  it  can  be  easily  removed  by  deconvolution. 

The  following  examples  demonstrate  the  result  of  using  the  above  two  $, 

approaches  for  removing  the  influence  of  the  driving  function  from  the 

response  function.  ►. 

Example  5 

Figure  3.13a  shows  the  response  function  starting  at  time  step  number 
sixty-one  and  taking  sixty  samples  at  every  third  time  step.  The  start- 
ing point  of  tg  » 61At  was  selected  for  two  reasons.  First  of  all,  by 
the  sixty-first  time  step  the  Gaussian  driving  function  value  is  well 
below  the  zero  value  of  the  computer.  Also,  as  was  mentioned  earlier, 
the  true  oscillatory  response  of  the  structure  does  not  begin  until  the 
current  has  propagated  to  the  end  of  the  structure  and  returned  to  the 
center  at  time  tg.  The  sampling  interval  of  3At  was  chosen  so  as  to 
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Figure  3.13.  (a)  Data  window  starting  at  time  step  sixty-one .Example  5. 

(b)  Resulting  pole  locations. 


49 


j ((i>L/C7T  ) 


I 


r ' 


i 


■ 

5 j 

Bi 

f 


El 


satisfy  the  Nyquist  criterion.  Since  sixty  samples  were  used,  the  method 
attempted  to  return  thirty  poles.  As  can  be  seen  from  Figure  3.13b,  only 
twenty-six  poles  were  actually  recovered,  four,  in  the  right  half  plane, 
were  thrown  away.  Of  the  twenty-six  returned,  only  twenty-two  are  the 
true  poles  of  the  system.  Thus,  even  though  the  driving  function  was 
removed,  extraneous  poles  are  still  present  due  to  the  fact  that  more 
than  the  twenty-two  true  poles  were  required  from  the  matrix  inversion 
process.  This  would  suggest  that  the  best  procedure  would  be  to  use  a 
least  squares  approach  and  search  for  only  twenty-two  poles. 

Example  6 

The  deconvolution  approach  is  demonstrated  in  Figure  3.14.  The  512 
samples  of  the  transient  response  have  been  transformed  to  the  spectral 
domain  using  a fast  Fourier  transform  routine.  This  spectral  response 
was  then  divided  by  the  spectral  response  of  the  Gaussian  pulse  and 
the  result  was  transformed  back  to  the  time  domain.  Figure  3.14a  shows 
the  first  eighty  samples  of  the  deconvolved  response  taken  at  a time  step 
size  of  2At.  Several  things  can  be  noticed  from  this  figure.  First  of 
all,  the  very  large  spike  at  about  t = 0 is  the  delta  function  which  is 
always  present  in  the  impulse  response.  Also,  high  frequency  ringing  is 
present  in  the  response  which  is  not  in  the  original  transient  response. 

This  high  frequency  ringing  is  due  to  Gibbs'  phenomenon  which  is  inherent 
in  the  deconvolution  process.  The  forty  poles  that  resulted  from  these 
eighty  time  samples  are  plotted  in  Figure  3.14b.  This  pole  pattern  is 
very  similar  to  that  of  Figure  3.4.  The  one  major  difference  is  apparent 
in  Figure  3.15,  which  is  a three-dimensional  plot  of  the  pole  pattern  show- 
ing the  amplitude.  Note  that  the  amplitude  of  the  true  poles  does  not  taper 
off  with  high  frequency  as  it  does  in  Figure  3.6  because  the  dependence 
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Figure  3.14.  (a)  Deconvolved  transient  response  data  of  Example  6. 

(b)  Locations  of  resulting  poles. 
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on  the  driving  function  has  been  removed  by  the  deconvolution  process. 

The  extraneous  poles  still  have  fairly  large  amplitudes  because  the  delta 
function  at  the  origin  was  included.  There  are  several  disadvantages 
to  performing  a deconvolution  to  remove  the  dependence  of  the  driving 
function.  The  main  disadvantage  is  the  large  amount  of  computation  time 
needed.  Several  stages  of  calculations  must  be  performed  which  cost 
time  and  accuracy  and  cause  errors  such  as  the  Gibb's  phenomenon.  In 
general,  the  deconvolution  process  should  be  avoided  if  at  all  possible. 

In  all  of  the  above  examples  the  extraneous  poles  were  present  be- 
cause the  number  of  time  samples  that  were  used  dictated  the  number  of 
poles  to  be  determined.  Also  in  the  above  examples,  errors  in  the  re- 
produced signals  resulted  because  in  calculating  the  N residues  only  N 
time  samples  were  used.  These  two  problems  can  be  overcome  by  using  the 
least-squares  error  approach  described  in  Section  3.1.  This  approach 
allows  one  to  use  many  time  samples  to  determine  just  a few  poles  thus 
permitting  one  to  select  beforehand  the  number  of  poles  to  be  determined. 
Examples  7 and  8 demonstrate  the  advantages  of  the  least-squares  or 
pseudo-inverse  solution  over  the  conventional  procedure. 

Example  7 

Figure  3.16a  is  a plot  of  the  sampling  interval  which  is  used  in 
this  example.  Note  that  120  samples  were  used  at  a sampling  rate  equal 
to  the  Nyquist  rate  for  the  problem,  3At.  The  first  twenty  samples  in 
the  interval  include  the  driving  function  so  the  resulting  poles  must 
be  expected  to  include  poles  representing  the  driving  function.  For  this 
example  only  twenty-six  poles  were  asked  for,  thus  allowing  for  the 
twenty-two  poles  which  were  known  to  be  the  true  poles  of  the  system  and 
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for  the  four  extra  poles  needed  to  represent  the  driving  function's  con- 
tribution. Figure  3.16b  is  a plot  of  these  poles  and,  indeed,  the  eleven 
true  pole  pairs,  which  were  expected,  appeared  along  with  two  extra  pole 
pairs.  The  interesting  point  here  is  that  although  only  twenty-six  poles 
and  residues  were  extracted  they  were  determined  using  a total  of  120 
samples  as  opposed  to  fifty-two  samples  which  would  have  been  used  if 
the  conventional  approach  had  been  employed.  Figure  3.17  demonstrates 
how  the  error  is  equally  distributed  over  the  entire  range  of  the  120 
samples  as  compared  to  the  error  of  Figure  3.10  which  is  near  zero  for  the 
first  twenty-eight  samples  and  then  increases  for  later  time. 

Example  8 

In  this  example,  Figure  3.18,  the  first  sixty  time  samples  have  been 
neglected.  Thus,  the  g(t)  term  is  eliminated  from  the  response  function 
and  only  the  true  poles  are  expected.  A total  of  150  samples  were  taken 
at  a rate  of  3At  and  only  twenty  poles  were  sought.  Figure  3.18a  shows 
that  the  resulting  poles  are  the  true  poles  of  the  system.  Note  from 
Figure  3.18a  that  as  frequency  increases  the  pole  pattern  diverges  from 
a sweeping  curve.  This  indicates  that  there  is  more  error  in  the  higher 
frequency  poles  due  to  the  lower  signal  level  at  those  frequencies. 

Note  that  the  frequency  u stays  stable  and  that  the  real  part  of  the 
pole  a is  most  sensitive  to  the  noise. 

Example  9 

In  all  of  the  previous  examples  transient  data,  which  were  generated 
using  a time-domain  computer,  were  used.  In  this  example  actual  experi- 
mental data  are  studied.  The  experimental  data  were  generated  on  the 
transient  electromagnetic  measurement  range  at  Lawrence  Livermore 
Laboratory  [11].  The  response  is  that  of  a 1.0  m monopole  located  on 
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a ground  plane  and  excited  by  an  approximation  to  a Gaussian-pulse  plane 
wave.  The  diameter  of  the  monopole  is  0.3175  cm.  The  antenna  is  loaded 
at  its  base  with  a 50  0 load  and  the  voltage  across  this  load  was  measured 
with  a sampling  oscilloscope.  A total  of  512  samples  were  measured  at  a 
time  interval  of  At  = 0.4  x 10  ^ s.  Of  the  512  measured  values  only  100 
samples  at  every  fifth-time  step  were  used.  Figure  3.19  shows  this 
measured  response  in  terms  of  the  current  through  the  load.  The  100  cur- 
rent samples  were  used  with  Prony's  method  and  forty-one  poles  and  residues 
were  produced.  Nine  poles,  which  were  in  the  right  half  plane,  were 
ignored.  Figure  3.20  is  a plot  of  the  generated  poles  found  in  the  second 
quadrant  of  the  complex  plane.  The  first  thing  that  is  apparent  about 
these  poles  is  that  they  fall  along  a curve  running  parallel  to  the  imag- 
inary axis.  This  is  typical  of  the  pole  locations  for  a dipole  as  seen 
from  the  previous  examples.  The  frequencies  of  the  first  nine  poles  in 
this  layer  correspond  to  the  first  nine  complex  resonant  frequencies  of 
aim  monopole.  The  real  parts  of  these  poles  seem  to  oscillate  around 
the  correct  value  because  of  the  sensitivity  of  the  real  part  to  the 
noise  in  the  data.  Even  though  the  response  was  measured  on  a sampling 
scope,  the  data  were  very  noisy  and  no  attempt  was  made  at  smoothing. 

The  fact  that  the  remaining  poles  do  not  correspond  to  physical  poles 
again  relates  to  the  fact  that  the  pulse  used  did  not  contain  frequency 
components  higher  than  that  of  the  ninth  resonance.  Also,  since  no 
attempt  was  made  to  remove  the  contribution  of  the  driving  function, 
these  extra  poles  are  needed  to  model  that  portion  of  the  response.  The 
pole  sitting  on  the  real  axis  close  to  the  origin  is  probably  due  to  the 
fact  that  there  was  a late  time  dc  level  present  due  to  the  pulser  used 
in  the  measurement  system. 
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3.4  Guidelines 


The  previous  section  presented  nine  examples  which  demonstrate  many 
of  the  aspects  of  Prony’s  method.  Contained  in  those  examples  are  guide- 
lines which  should  be  carefully  followed  when  using  the  method.  Since 
the  guidelines  were  somewhat  obscured  by  the  numerous  examples,  they 
will  be  repeated  in  this  section  as  a summary  for  ready  reference. 

3.4.1  Sampling  interval 

The  width  of  the  sampling  interval  should  be  selected  so  that  it 
is  wide  enough  to  contain  all  of  the  characteristics  of  the  response. 

For  the  most  part,  the  interval  should  be  at  least  as  long  as  one  period 
of  the  lowest  frequency  present.  If  a least  squares  method  is  employed, 
then  the  length  of  the  sampling  interval  should  be  as  long  as  is  eco- 
nomically feasible. 

3.4.2  Sampling  rate 

All  physically  obtainable  transient  responses  will  be  bandlimited 
and  thus  the  upper  frequency  limit  should  be  determinable.  The  sampling 
rate  should  then  be,  by  the  Nyquist  criterion,  slightly  more  than  twice 
the  highest  frequency  expected. 

3.4.3  Removal  of  the  influence  of  the  driving  function 

The  influence  of  the  driving  function  on  the  response  function  should 
be  removed  or  accounted  for  by  using  one  of  the  processes  described  in 
Section  2.2.  The  two  easiest  processes  to  use  consist  of  either  a driv- 
ing function  that  can  be  represented  by  known  poles  or  a driving 
function  which  is  time  limited.  If  a time-limited  driving  function  is 
used,  sufficient  response  samples  must  occur  after  the  turnoff  time  of 
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the  driving  function  so  that  the  required  number  of  poles  may  be  determined. 

In  an  experimental  system  where  this  may  not  be  practical,  the  driving 
function  should  be  representable  by  a sum  of  sinusoids.  Deconvolution 
should  be  avoided  because  of  the  expense  and  the  increase  in  the  error 
level. 

3.4.4  Least  squares  versus  conventional  matrix  inversion 

* 

The  least-squared  error  or  pseudo- inverse  approach  should  always  be 
used.  This  allows  as  many  samples  to  be  used  as  one  desires  and  allows 
for  the  selection  of  the  number  of  poles  without  being  dependent  on  the 
number  of  samples  used.  This  method  gives  better  results  in  calculating 
the  residues  since  more  input  samples  can  be  used. 

3.5  Problems  Associated  with  Prony's  Method 

There  are  two  major  problems  associated  with  Prony's  method  which 
need  to  be  overcome  in  order  for  the  method  to  be  practical.  The  first 

problem  is  the  fact  that  it  is  necessary  to  know  the  number  of  poles  N ' 

which  are  contained  in  the  transient  data  before  Prony's  algorithm  can 
be  applied.  The  second  problem  is  that  Prony's  method  is  extremely 
sensitive  to  noise  of  any  kind.  These  problems  are  discussed  here  and 
solutions  to  these  problems  are  presented  in  the  next  two  chapters. 

The  numerical  examples  of  Section  3.3  indicate  that  if  one  asks 
for  more  poles  than  are  effectively  contained  in  the  data  then  the 
algorithm  generates  a number  of  extraneous  poles  in  addition  to  the  ones 
that  compare  to  the  true  poles.  The  presence  of  the  extraneous  poles 
causes  the  residues  of  the  true  poles  to  be  inaccurate  and  results  in 
unnecessary  computation  time.  Similarly,  if  one  underestimates  the 
number  of  poles, then  many  of  the  returned  poles  may  substantially 
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deviate  from  the  true  poles  and,  in  most  cases,  the  true  poles  will  not 


be  returned  at  all.  Thus,  a systematic  approach  for  an  a priori  deter- 
mination of  N would  be  extremely  beneficial.  The  examples  showed  that, 
if  the  upper  frequency  limit  is  known  and  if  the  approximate  values  of 
the  poles  are  known,  one  can  make  a good  estimate  of  the  value  of  N. 

Also  N could  be  determined  by  trial  and  error  but  that  would  be  very 
expensive.  The  next  chapter  presents  two  straightforward  and  systematic 
schemes  for  determining  N. 

The  other  problem  with  Prony's  method  is  the  sensitivity  of  the 
poles  to  a noisy  set  of  data.  The  real  part  of  the  poles  is  extremely 
sensitive  to  noise.  If  the  noise  is  too  bad,  then  Prony's  method  will 
not  return  any  of  the  true  poles  and  will  attempt  to  curve-fit  to  the 
noise.  Therefore,  it  is  necessary  to  determine  the  noise  limitations 
for  Prony's  algorithm  and  to  determine  ways  of  preventing  and  reducing 
noise.  This  problem  is  handled  extensively  in  Chapter  5. 
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4.  THE  DETERMINATION  OF  THE  NUMBER  OF  POLES  (N) 


In  the  last  chapter  it  became  obvious  that  in  order  to  use  Prony's 
method  optimally  a systematic  approach  must  be  developed  for  determining 
the  number  of  poles  contained  in  the  transient  response.  This  chapter 
develops  two  techniques  by  which  the  number  of  poles  inherent  in  the 
system  can  be  analytically  determined.  The  first  approach  discussed 
was  developed  by  Householder  [27].  His  derivation  is  presented  here  ex- 
cept that  it  is  applied  to  the  least  squares  solution.  The  second  ap- 
proach is  new  and  is  shown  to  have  several  advantages  over  Householder's 
method. 


4.1  Householder  Orthogonalization  Method 

It  was  shown  in  Chapter  3 that  if  a set  of  discrete  transient  re- 
sponse data  could  be  represented  as  a sum  of  N complex  exponentials  then 
the  discrete  samples  1^  must  satisfy  a difference  equation  of  order  N. 
Consider  now  the  following  vectors  which  are  filled  with  the  data 
samples  I. 


X0  ’ l2  ’ X1  ” *3 


(4.1) 


where  the  ith  vector  is  the  i - 1th  column  of  matrix  (3.8b)  in  Chapter  3. 
Since  the  1^  satisfy  Prony's  difference  equation  of  order  N,  (3.5b), 
then  so  must  the  vectors  i^.  Now  if  y > N,  then  any  N,  or  fewer,  of  the 
vectors  i^  will  be  linearly  independent  and  any  N + 1 of  the  vectors  will 
be  linearly  dependent.  Hence,  the  difference  equation  can  be  written 


i0  “0  + il  al  + • 


+ ^ aN 


(4.2) 


The  next  step  is  to  apply  the  Gram  Schmidt  or  Choleski  orthogonalization 
process  to  the  vectors  i^. 

The  orthogonalization  process  consists  of  replacing  each  vector  i^ 
by  a linear  combination  of  the  vectors  iQ,i^,  • • •,  i^.  The  new  ortho- 
gonal vectors  are  denoted  as  a^a^,  • • . , a^  The  steps  are  then 

a0  “ x0 

ai  * h - vio  ao 

where  is  chosen  so  as  to  make  a^  and  a^  orthogonal.  That  is, 

T 


10 


iiao 


aoao 


where  an  is  the  transpose  of  vector  an.  Then  proceeding  to  the  second 


vector 


a2  * *2  " U21  al  ' g20  a0 


where  again  and  i^q  are  selected  so  that  a2  is  orthogonal  to  both  a^. 


and  a^,  which  means  that 
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a0  i *2 

J20  = ~T~ 
i i 
0 10 


T. 

, - 

21  T 

i f 

xi  h 


Continuing  this  process,  one  finally  obtains 


Vl  ■ Vl  " PN-l,N-2  aN-2  + 


* + UN-1,0  30 


where 


^-2  Sj-l 
N-l.N-2  “ T 

aN-2  aN-2 


ao  Vl 


N-1,0  T 


ao  a0 


Let  V1,  and  designate  the  matrices 


Vi  = ‘ * - * V 


Aj  - ^a0,al’  . . a^) 


1 U10  u20  * 


0 1 


(4.5) 
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Then, 


V ■ A M 
i i i 


(4.6) 


which  is  the  Gram-Schmidt  procedure  in  matrix  form.  Since  is  a 
triangular  matrix,  it  is  easily  inverted  so  that 


A.  - V,  M. 
i i l 


The  difference  Equation  (4.2)  can  now  be  written  as 


VN-1  °N  “ _1N 


(4.9) 


Equation  (4.9)  is  equivalent  to  Equation  (3.6)  of  Chapter  3.  Now  if  a 
pseudo-inverse  approach  is  used  to  solve  Equation  (4.8)  then  the  re- 
sulting expression  is 


v»  i . o„  - vl  , (-i  ) 
N-l  N-l  N N-l  n 


(4.10) 


From  (4.6) 


VN-1  ” 'Vl  ^J-l 


(4.11a) 


VN-1  * *4-1  Vi 


(4.11b) 
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so  that  (4.10)  can  be  written  as 


^-1  ^-1  "Sj-l  ^-1  °N  “ ^-1  V-l^V 


(4.12) 


\-l  'Vl  = °N-1 


(4.13) 


T 

where  the  matrix  is  a simple  diagonal  matrix  because  A^  ^ and  A^ 

are  made  up  of  orthogonal  vectors.  Now  define  the  vector 


Ji  - ‘ 


(4.14) 


Ji,i-1 


so  that  the  orthogonalization  process  can  be  written 


Li+1  = ai+l  + Ai  Ui+1 


(4.15) 


or  by  multiplying  both  sides  by  A* 


Ai  ii+l  “ Ai  ai+l  + Ai  Ai  ui+l 


(4.16) 


Since  all  columns  of  A^  are  orthogonal  to  a^+^,  then 


Ai  ii+l  ” Di  wi+l 


(4.17) 
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where  has  replaced  the  product  A^.  If  Equations  (4.17)  and  (4.13) 
are  substituted  in  (4.12),  then 


“n-I  °N 


UN(-i) 


(4.18) 


Since  the  bracketed  terms  on  each  side  of  (4.18)  are  nonsingular  matrices, 
they  can  be  removed  by  multiplying  both  sides  by  their  inverse.  Hence, 
(4.18)  reduces  to 


^-1  °N  ’ -UN  ' (4,1 

To  summarize,  the  and  a^+^  are  calculated  sequentially  by  using 

the  Gram-Schmidt  procedure.  Equations  (4.15) and  (4.17).  The  vector 

a^+^  is  adjoined  to  A^  to  obtain  A^+^  and  matrix  is  bordered  by 

and  a unit  row  vector  to  obtain  M.,,.  Note  that  each  vector  aJ  . is 

l+l  i+1 

the  component  of  i^+^  that  is  orthogonal  to  the  space  of  the  previous 
i's.  Since  the  vectors  i^  must  satisfy  a difference  equation  of  order 
N,  then  i^  is  a linear  combination  of  the  preceding  i's  and  a^  should 
vanish.  Thus,  in  order  to  determine  the  value  of  N,  the  above  process 
should  be  continued  until  one  of  the  vectors  a vanishes.  However, 
since  the  time  samples  are  subject  to  measurement  errors  and  roundoff 
errors,  it  should  not  be  expected  that  any  vector  a will  vanish  entirely 
but  for  some  N there  will  be  an  which  will  be  negligibly  small. 

Once  a negligibly  small  vector  a^  has  been  found,  it  is  a simple 
procedure  to  find  the  coefficients  a of  the  difference  equation.  The 
coefficients  a are  found  from  (4.19).  Note  that  the  matrix  is  a 

triangular  matrix  so  that  the  vector  may  be  found  by  back  substitution, 
eliminating  a matrix  inversion.  Thus,  this  process  not  only  gives  a 
method  for  determining  the  number  of  poles  N in  the  system  but  also 


K 


allows  one  to  actually  find  the  N poles  without  having  to  invert  a 
matrix.  Note  also  that  this  method  can  use  as  many  time  samples  as  de- 
sired by  selecting  the  appropriate  value  of  y in  (4.1). 

As  an  example  of  this  procedure,  a set  of  test  transient  data  was 
produced  using  the  set  of  six  pole  pairs  listed  in  Table  4.1.  For  all 
the  poles  the  accompanying  residues  were  taken  to  be  one.  The  resulting 
transient  response  is  plotted  in  Figure  4.1.  Householder's  method  was 
carried  out  exactly  as  outlined  above.  The  value  of  y which  was  used 
was  100.  Figure  4.2a  is  a plot  of  the  average  of  the  absolute  value 
of  the  first  thirteen  vectors  i^  and  the  resulting  orthogonal  vectors 
a^.  Note  that  while  the  average  of  the  first  thirteen  vectors  i^  stays 
constant  the  first  twelve  values  of  a^  drop  off  slightly  until  at  the 
thirteenth  vector  the  value  drops  by  four  orders  of  magnitude.  While 

this  thirteenth  vector  did  not  go  to  zero,  it  did  drop  low  enough  to 

< 

indicate  that  the  thirteenth  vector  was  linearly  dependent  on  the  pre- 
ceding twelve  vectors.  The  difference  equation  coefficients  were  calcu- 
lated as  per  Equation  (4.19)  and  the  resulting  poles  are  listed  in 
Column  Two  of  Table  4.1.  Note  that  the  resulting  poles  are  essentially 
the  same  as  the  original. 

As  another  test  of  the  above  procedure,  the  same  transient  response 
was  used  but  normally  distributed  noise  with  a standard  deviation  of 
0.005  was  added  to  the  response.  Again  y was  set  equal  to  100  and  the 
resulting  vectors  a^  and  i^  are  plotted  in  Figure  4.2b.  Notice  now 
that  at  the  thirteenth  vector  the  average  value  does  not  drop  off.  This 
indicates  that  the  Householder  procedure  is  not  applicable  to  systems 
with  significant  noise.  When  the  poles  were  found  for  the  noise  case, 
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TABLE  4 


NUMBER  OF  VECTOR 


NUMBER  OF  VECTOR 

(b) 

Figure  4.2.  Value  of  the  average  of  the  absolute  value 
of  the  vectors  a and  i. 

(a)  No-noise  case. 

(b)  Noise  of  0 = 0.005  added  to  data. 
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they  were  not  associated  with  the  true  poles  at  all,  as  indicated  by  the 
third  column  of  Table  4.1.  Chapter  5 shows  that  the  alternate  procedure 
for  finding  the  poles,  which  is  presented  in  the  next  section,  does  in- 
deed work  in  the  presence  of  noise. 

As  a further  example  the  numerically  generated  transient  response 
of  the  current  on  the  one  meter  dipole  discussed  in  Section  3.3  is  used. 
The  samples  were  taken  at  every  third  time  step  starting  with  time  step 
sixty-one.  Figure  4.3a  shows  a plot  of  both  vectors  i and  the  resulting 
orthogonal  vector  a.  Note  that,  while  vector  i has  an  average  which  is 
somewhat  constant,  vector  a drops  off  until  vector  number  twenty-three 
where  it  stays  level  for  a while.  This  is  interesting  since  in  the  ex- 
amples of  Section  3.3  it  was  shown  that  eleven  pole  pairs  appeared  to 
be  the  optimal  solution  to  this  example.  Indeed,  when  the  process  of 
(4.19)  was  applied  using  the  twenty-third  orthogonal  vector,  the  first 
eleven  poles  in  Figure  3.4  were  obtained.  Note  that  the  twenty-third 
vector  a is  four  orders  of  magnitude  below  the  twenty-third  vector  i. 

In  Figure  4.2a  of  the  previous  example, the  thirteenth  vector  a was  also 
four  orders  of  magnitude  below  the  thirteenth  vector  i.  This  implies 
that  when  the  average  of  the  orthogonal  vector  drops  four  orders  of 
magnitude  below  its  accompanying  response  vector  then  that  vector  is 
the  first  dependent  vector. 

In  Figure  4. 3b  the  vectors  a and  i are  plotted  for  the  dipole  case 
.nsldered  above, but  since  the  samples  start  at  time  step  one,  the 
:n<  function  is  included.  Note  that  none  of  the  twenty-nine  a 

ire  -nore  than  one  order  of  magnitude  below  the  vector  i.  This 
•r«  -i  <:nce  the  Gaussian  pulse  included  in  the  first  sixty  samples 
• - 4*  a finite  sum  of  exponentials. 
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Figure  4.3.  Value  of  Che  average  of  the  absolute  value  of  the 

vector  a and  i for  the  response  of  the  1.0  m dipole 

(a)  Data  window  starting  at  tine  step  sixty-one. 

(b)  Data  window  starting  at  tine  step  one. 


To  start  this  development,  it  is  again  necessary  to  form  the  vectors 


vi  of  Section  4.1.  Note  that  the  v^  can  be  written  as 
N s nAt 

v = £ A e ^ n - i,i+l,  . . i + y 

J-l  J 

or  by  separating  terms 

N s iAt  s nAt 

v,  = I (A4  e ^ ) e J n = 0,1,  ... 

J-l  J 

Expression  (4.21)  can  be  written  as 


v . 
i 


N 

■ l 

J-l 


“J.i  "J 


(4. 


(4.; 


(4. 


where 


C3,i  ‘ Ai  * 


Sj  iAt 


(4. 


t h th 

is  the  coefficient  for  the  j mode  vector  i(i  starting  at  the  i time 


step.  Remembering  that 


v 


Sj  At 


,th 


then  the  i mode  vector  is 


(4. 


If  there  are  N poles  in  the  system,  then  there  will  be  N mode  vectors 


The  next  step  is  to  solve  the  difference  equation 


N 

l V ..  a 
m-0  m+k  ra 


0 


(4.26) 


If  a pseudo- inverse  solution  is  used,  then  (4.26)  would  appear  in  matrix 
form  as 


(4.26) 


Matrix  $ is  defined  as  the  product  of  the  above  two  matrices.  Note  that 
the  it'1,  jCk  element  of  the  $ matrix  is 


i.J 


(4.27) 


T 

where  v^  is  the  transpose  of  vector  v^  Substituting  (4.23)  into  (4.27) 
gives 


N N 


- I l 

£**1  m«l 


T* 

\ *m 


(4.28) 


where  the  * indicates  complex  conjugate.  Since  $ is  a matrix  of  order 
S + 1 by  K f 1,  it  will  have  one  eigenvector  which  will  be  orthogonal 
to  the  N mode  vectors.  That  is,  there  will  be  one  eigenvector  such  that 


4>  E 


N+l 


0 


(4.29) 
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where  is  the  eigenvector  orthogonal  to  all  of  the  modal  vectors. 

The  eigenvalue  corresponding  to  this  eigenvector  is  zero.  If  matrix 
# was  made  to  be  2N  there  would  be  N eigenvectors  orthogonal  to  the  N 
modal  vectors  and  there  would  be  N eigenvalues  equal  to  zero.  Hence, 
the  process  for  determining  the  N of  the  system  is  to  fill  matrix  $ 
to  some  dimension  M by  M.  The  corresponding  M eigenvalues  of  the 
system  are  found  and  checked  to  see  if  one  or  more  is  equal  to  zero. 

If  there  are  L eigenvalues  equal  to  zero,  then  N would  be  equal  to  M - L. 
If  L is  not  equal  to  one,  then  the  matrix  $ is  recomputed  to  order 
N + 1 by  N + 1 and  the  eigenvalues  regenerated.  Equations  (4.26)  and 
(4.29)  show  that  the  eigenvector  corresponding  to  the  one  zero  eigen- 
value is  the  vector  of  the  coefficients  of  the  difference  equation. 

Thus,  not  only  are  the  number  of  poles  found,  but  the  poles  are  found 
at  the  same  time. 

As  the  first  example  of  this  procedure,  consider  again  the  tran- 
sient response  of  Figure  4.1  which  was  generated  using  the  six  pole 
pairs  of  Column  One  of  Table  4.2.  The  thirteen  eigenvalues  of  this 
system,  which  were  generated  using  a value  of  y = 100,  are  plotted  in 
Figure  4.4.  Notice  the  very  sharp  drop  between  the  number  twelve  and 
the  number  thirteen  eigenvalue.  The  twelve  poles  that  resulted  from 
the  eigenvector  corresponding  to  the  thirteenth  eigenvalue  are  listed 
in  Column  Two  of  Table  4.2.  The  resulting  poles  are  in  good  agreement 
with  the  original. 

As  a further  test  of  the  above  procedure,  the  same  transient  data 
were  used  but  normally  distributed  noise  with  a standard  deviation  of 
0.005  was  added.  The  value  of  y was  set  to  100  and  the  resulting 
thirteen  eigenvalues  are  plotted  in  Figure  4.4.  Note  that  the  drop 
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between  the  twelfth  and  thirteenth  eigenvalues  is  only  slightly  more  than 
one  order  of  magnitude,  indicating  that  the  noise  has  a gross  effect  on 
this  procedure.  The  poles  that  were  recovered  from  this  noisy  case  are 
listed  in  Cr’.umn  Three  of  Table  4.2.  Although  the  values  of  these  poles 
are  not  precisely  the  true  poles,  the  maximum  percent  error  is  only  1.5. 
When  this  fact  is  compared  to  the  fact  that  some  of  the  poles  of  Column 
Three  of  Table  4.1  are  not  even  in  complex  conjugate  pairs,  the  conclusion 
indicates  that  this  procedure  will  work  with  noisy  data.  The  next  chapter 
will  show  that  the  eigenvalues  are  also  a function  of  the  noise 
level. 

As  a final  example  of  this  procedure  again  consider  the  transient 
response  data  of  Figure  4.1.  Samples  were  taken  at  every  third  time 
step  starting  at  time  step  sixty-one.  Tests  were  run  in  which  twenty- 
one,  twenty-three,  and  twenty-five  eigenvalues  were  determined  and  the 
corresponding  poles  were  extracted.  Figure  4.5  shows  plots  of  the 
twenty-three  eigenvalues  and  twenty-five  eigenvalues  and  Figure  4.6 
plots  the  poles  for  the  twenty,  twenty-two,  and  twenty-four  pole  cases. 

For  the  case  of  twenty-two  poles,  the  poles  correspond  to  the  first 
eleven  poles  plotted  in  Figure  3.4,  while  for  the  case  of  twenty  poles 
the  upper  frequency  poles  vary  significantly.  The  case  where  twenty- 
four  poles  were  evaluated  gives  two  poles  that  do  not  form  complex  con- 
jugate pairs  and  the  eleventh  pole  pair  does  not  correspond  to  the 
eleventh  pole  pair  of  the  twenty-two  pole  example.  This  indicates  that 
the  case  where  twenty-three  eigenvalues  were  generated  is  the  proper 
solution  and  should  be  studied.  Notice  that  in  the  plot  of  the  twenty- 
three  eigenvalues  the  eigenvalues  drop  off  continuously  until  the 
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Figure  4.5.  Comparison  of  the  resulting  twenty-three  and 

twenty-five  eigenvalues  for  1.0  m dipole  response. 
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twenty-third  without  any  significant  jump,  except  possibly  between 
twenty-two  and  twenty-three.  This  is  probably  due  to  the  fact  that  the 
residues  of  these  eleven  pole  pairs  drop  off  significantly  as  frequency 
increases.  Because  of  this  lack  of  a significant  breakpoint,  it  appears 
that  this  technique  is  not  as  well  suited  to  determining  the  value  of  N 
for  systems  which  have  transient  responses  similar  to  the  response  of 
this  example. 
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5.  NOISE  AND  ITS  RELATIONSHIP  TO  PRONY'S  METHOD 

In  the  preceding  chapters  several  allusions  were  made  to  the  fact 
that  noise  seriously  affects  the  values  of  the  poles  extracted  by  Prony's 
method.  In  Chapter  3 it  was  pointed  out  that  if  the  noise  level  is  high 
enough  Prony's  method  will  return  poles  that  are  not  even  remotely  as- 
sociated with  the  true  poles  being  sought.  This  effect  can  be  understood 
if  the  noise  is  thought  of  as  being  several  arbitrary  frequency  com- 
ponents added  to  the  signal.  Thus,  the  true  exponential  nature  of  the 
signal  has  been  corrupted.  Since  Prony's  method  is  an  interpolation  pro- 
cess, it  will  give  a set  of  poles  which  fit  the  noisy  transient  response 
but  will  not  necessarily  be  the  natural  resonances.  The  least  squares 
approach  is  used  to  reduce  some  of  the  error  and  to  allow  more  data 
samples  to  be  used.  This  chapter  shows  that  the  least  squares  approach 
when  applied  to  Prony's  method  does  not  strictly  give  a least  squares 
fit.  If  the  noise  level  is  known,  it  can  be  used  to  aid  in  the  deter- 
mination of  the  number  of  poles  N using  the  eigenvalue  process  of 
Chapter  4.  Statistical  studies  are  presented  which  relate  the  noise 
levels  to  the  quality  of  the  results  obtained. 

5. 1 Least  Mean  Squared  Error  Approach 

Up  to  this  point  the  least  squares  approach  has  been  applied  blindly 
by  simply  performing  a pseudo-inverse  solution.  This  was  done  primarily 
to  allow  for  more  transient  data  points  to  be  used  when  solving  for  a 
set  number  of  poles.  This  section  studies  the  way  in  which  the  least 
mean  squared  error  process  reduces  the  errors  that  are  present  in  the 


transient  data. 


Consider  a set  of  transient  response  data  which  has  been  sampled 


at  M equally  spaced  intervals  At.  These  M sampled  values  Rg,  R1 , 
are  assumed  to  be  estimates  of  the  true  response  of  the  system 


R0’  Rl’  * 


where 


s ,nAt 


R = R(nAt)  = £ A.  e * n = 0,  1,  . . M 


(5.1) 


The  measured  sampled  values  R differ  from  the  true  values  R because 

n n 

of  errors  in  measurements,  noise,  etc.  It  is  known  that  the  R satisfy 

I\ 

Prony's  difference  equation  of  order  N which  may  be  written 

Rn+N  + aN-l  Rn+N-1  + * * ’ + ai  Rn+1  + a0  Rn  = °;  n = °»1’  * ' * ’ M ' 

(5.2) 

In  Chapter  3 it  was  shown  that  the  pseudo- inverse  or  least  squares  solution 
to  (5.2)  could  be  written  in  matrix  form  as  (3.8a) 

AT  A B = AT  C (5.3) 


4>  B = D 


(5.4a) 


where 


l Rj  l Rj  R. 

j=*0  J j=0  J • 


■ Jo  Rj  Vn-1 


$ * A A 


Rj+1  Rj 


l . 

* ’jio  Vl  Rj+N-1 


X . I . . 1 . 

Rj+n-i  Rj  4L  rj+n*1  Vi  ’ ' 


(5.4b) 


The  y of  (5.4b)  and  (5.4d)  must  be  greater  than  2N  - 1 and  less  than 

M - N - 1.  Even  though  Equations  (5.4)  are  formed  in  accordance  with 

the  usual  least  squares  approach,  the  usual  assumptions  do  not  hold. 

Since  the  are  subject  to  errors,  both  the  left  side  and  the  right 

side  of  (5.4a)  are  formed  by  noisy  elements.  In  the  normal  least 

squares  approach,  only  the  unknown  quantities  are  assumed  subject  to 

error.  Thus,  when  the  pseudo-inverse  procedure  is  applied  to  Prony's 

method,  it  does  not  yield  an  optimum  least  squares  approximation  since 

the  M samples  R do  not  represent  exactly  a sum  of  N exponentials; 
n 

that  is, 


(5.4c) 


(5.4d) 
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where  is  the  error  in  the  n measured  time  sample.  If  equals 
zero,  then  the  pseudo-inverse  is  the  optimum  least  squares  approximation. 

Householder  [27]  developed  a method  by  which  the  N parameters 
could  be  determined  such  that  the  sum  of  the  squared  error 


S = J w (R  - R ) 
L n n " 


(5.6) 


is  minimized  subject  to  the  fulfillment  of  the  side  condition  (5.2). 

The  multipliers  w^  are  the  statistical  weights  associated  with  the 
measured  values  R^.  Householder's  procedure  is  an  iterative  procedure 
involving  Lagrange  multipliers.  The  method  is  repeated  several  times 
until  the  results  are  within  the  accepted  tolerance  level.  The  method 
is  very  laborious,  and  if  more  than  just  a few  poles  are  sought,  it  is 
totally  unwieldy.  In  1966,  McBride,  Schaefgen,  and  Steighlitz  [28]  and, 
in  1968,  McDonough  and  Huggins  [16]  developed  two  different  linear  iter- 
ative schemes  which  appear  quite  successful  for  determining  the  poles 
even  for  large  N.  The  methods,  however,  were  developed  for  synthesizing 
a prescribed  transient  response  using  a sum  of  N complex  exponentials. 

It  was  not  required  that  the  poles  found  be  the  true  waveform  poles, or 
indeed,  it  was  not  required  that  the  waveform  even  have  N poles.  Hence, 
these  methods  are  also  not  satisfactory  for  the  requirements  here,  that 
is,  finding  the  N true  poles  from  a set  of  noisy  transient  response  data 
of  a system  having  precisely  N poles. 

The  next  section  presents  a method  by  which  the  poles  of  a system 
can  be  extracted  from  noisy  data  using  an  eigenvalue  approach  similar 
to  that  presented  in  Section  4.2. 


- 
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5.2  Determination  of  the  Poles  from  Noisy  Data 


In  this  section  the  approach  in  Section  4.2  is  applied  to  tran- 
sient responses  containing  noisy  data.  Consider  the  transient  response 
vector  of  (4.23)  to  include  noise  such  that 


'i  -Ji  cj.i  *i + ni 


(5.7) 


where  is  the  noise  vector  starting  at  the  iCh  time  step. 


'i+1 


(5.8) 


"Y+i 


Assume  that  the  values  of  noise  vector  are  due  to  a random  process 

2 

and  are  normally  distributed  with  zero  mean  and  variance  a . The 

definitions  of  iji  , the  jtt;  mode  vector,  and  the  coefficients  C . are 
J J » 1 

precisely  those  of  (4.25)  and  (4.24),  respectively.  As  before,  if 
there  are  N poles  in  the  system  then  there  will  be  N mode  vectors  ifi.  . 

The  next  step,  as  in  Section  4.2,  is  to  solve  the  difference  equation 


N 


£ v .,_  a_  ■ 0 


m-»0 


nt+k  m 


(5.9) 


If  a pseudo-inverse  solution  of  (5.9)  is  attempted,  the  set  of  Equations 
(4.26)  is  obtained.  If  the  noisy  transient  response  vectors  are 


I 


substituted  into  the  matrix  4>  of  (4.27),  the  following  result  is  obtained 


for  the  i,j  element  of  matrix  <t>',  where  the  prime  indicates  that  the 


matrix  was  formed  using  the  noise  vector 


t>'  = VT  v 

i.j  i J 


(5.10a) 


N N 


*ll1  [c>-.j  ni'  *i + lni  "j1  • 


(5.10b) 


Since  the  term  [r^  n j ] is  the  product  of  the  transpose  of  the  ith  noise 


vector  times  the  j noise  vector  and  since  any  two  noise  vectors  are 


uncorrelated,  then 


£n^  Hj]  = 0 , i 4 j 


(5.11a) 


= Y a , i = j 


(5.11b) 


where  y is  the  length  of  the  response  vector  v^  and  the  length  of  the 


noise  vector  • Hence,  this  term  gives  rise  to  a matrix  which  is  zero 


everywhere  and  has  a value  of  y o on  the  diagonal,  or  more  precisely 


written. 


H - y a I 


(5.12) 


where  I represents  the  identity  matrix.  The  second  and  third  terms  on 


the  right  side  of  (5.10b)  are  each  zero  since  they  are  the  product  of 


the  noise  vector  and  the  natural  mode  vector  which  are  uncorrelated  and 


r 


of  zero  mean.  The  first  term  is  left  which  is  precisely  the  <f  of  (4.28). 

th 


This  term  is  the  product  of  the  l and  n natural  mode  vectors.  Thus, 

- $ + H (5.13) 

where  <t> ' is  of  order  N + 1 by  N + 1.  Matrix  4>'  will  have  N + l real 
eigenvectors  of  which  one  eigenvector  will  be  orthogonal  to  the  N mode 


vectors  . That  is, 
m 


Et  * 0 , l - l.N 


0 , l > N 


(5.14) 


where  E,  is  the  transpose  of  the  ith  eigenvector.  Hence, 
E. 


“N+l  * Y 0 EN+1 


(5.15) 


since  4>  E^+^  ■ 0 and  matrix  H can  be  thought  of  as  the  constant  y o . 

, th 


Thus,  (5.15)  implies  that  the  eigenvalue  associated  with  the  N+l 


eigenvector  is  Just  y •;  which  is  y times  the  variance  of  the  noise. 


This  is  similar  to  the  result  of  Section  4.2  where  the  eigenvalue 

.th 


associated  with  the  N+l  eigenvector  is  zero  for  the  noise-free  case. 
The  procedure  for  determining  the  value  of  N is  then  the  same  as  that 


outlined  in  Section  4.2  except  that  matrix  is  filled  until  it  has  an 


2 th 

eigenvalue  equal  to  y o instead  of  zero.  Likewise  the  N+l  eigen- 


value, which  is  orthogonal  to  all  the  mode  vectors,  is  the  vector  con- 


taining the  N+l  coefficients  of  Prony's  difference  equation.  Thus, 

2 


once  the  eigenvalue  equal  to  y o is  found,  the  poles  can  be  found  from 
its  corresponding  eigenvector. 


, 
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The  next  section  applies  this  procedure  to  several  examples. 
Different  noise  levels  are  used  in  order  to  statistically  determine  the 


L 


l 


poles'  sensitivities  to  various  noise  variances. 

5.3  Numerical  Examples 

Example  1 

As  the  first  example  of  the  method  developed  in  Section  5.2,  a 
single  undamped  sinusoid  is  used.  The  transient  data  were  produced  at 
200  time  steps  using  the  relation 

R(nAt)  = sin(TrnAt);  n = 0,1,  . . .,  199 

where  At  » 0.1  s.  Noise  of  different  levels  was  added  to  the  data  samples 
R(nAt).  The  noise  was  produced  by  using  a pseudo-random  number  gen'  ator 
on  a digital  computer  to  produce  uniformly  distributed  noise.  The  uni- 
formly distributed  noise  was  then  transformed  into  normally  distributed 
noise  with  zero  mean  and  a standard  deviation  of  o.  For  each  different 
value  of  a,  the  entire  pole  extracting  process  was  repeated  twenty  times 
with  twenty  different  sets  of  random  numbers.  That  is,  twenty  Monte  Carlo 
trials  were  performed  for  each  standard  deviation  of  the  noise.  From 
these  twenty  trials  expected  values  of  the  poles  were  calculated  along 
with  the  variance  of  the  poles.  The  expected  value  of  the  N + ] eigen- 
value and  its  variance  were  also  calculated.  Table  5.1  shows  the  percent 
error  of  both  the  real  and  imaginary  parts  of  the  poles  produced  for 
different  values  of  a.  The  percent  error  was  calculated  as 


M 

I 


[1 
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where  ST  is  the  true  value  of  the  complex  pole  (in  this  case  ST  » 0 + jtr) 
and  SE  is  the  calculated  expected  value  of  the  extracted  poles  for  a 
particular  noise  level.  Also  tabulated  in  Table  5.1  are  the  standard 
deviation  of  the  extracted  poles,  the  theoretical  value  of  the  N + 1th 
eigenvalue  (i.e.,  y o^) , the  mean  value  of  the  calculated  N + 1C^  eigen- 
value, and  the  signal-to-noise  level  in  decibels.  The  signal-to-noise 
level  was  calculated  by  taking  the  log  of  the  square  of  the  average  of 

the  absolute  value  of  the  signal  over  the  time  window  of  length  yAt 

2 

and  dividing  by  the  variance  of  the  noise  o . That  is 

( Y-  I R(iAt)  1 

S/N (dB)  - 10  log  

o 

For  this  first  example  the  average  of  the  absolute  value  of  the  signal 
is  0.63137  for  all  time  because  the  signal  is  undamped.  The  values  of  a 
for  the  noise  used  ranged  from  0.5  to  0.01  or  in  terms  of  signal-to-noise 
level  from  2.0  dB  to  36.0  dB.  Note  that  the  percent  error  of  the  real 
part  of  the  poles  is  more  sensitive  to  noise  than  is  the  imaginary  part. 
Note  also  that  the  theoretically  predicted  value  of  the  N + 1^  eigen- 
value is  extremely  close  to  the  calculated  value  of  that  eigenvalue. 

This  example  indicates  that  the  poles  can  be  extracted  from  data  when 
the  signal-to-noise  level  is  as  low  as  2.0  dB.  If  the  standard  deviation 
of  the  extracted  poles  is  studied  for  the  2.0  dB  case,  the  results  do 
not  look  very  good.  Note  that  for  the  imaginary  part  of  the  pole  the 
standard  deviation  is  0.796  which  is  extremely  high  when  it  is  remembered 
the  value  of  the  pole  is  just  n.  The  standard  deviation  of  the  pole  at 
the  16.0  dB  level,  however,  is  much  smaller  indicating  that  the  chance  of 


. 1 


I 

t 
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extracting  the  true  pole  is  much  higher  at  16.0  dB  than  at  2.0  dB.  This 
result  will  be  further  substantiated  in  the  next  example. 

Example  2 

For  the  second  example  the  same  sinusoid  as  in  Example  1 was  used 
but  now  it  is  exponentially  damped.  That  is, 

R(nAt)  = e 2,0nAt:  sin(nnAt)  n = 0,1,  . . . , 199 

where  again  At  = 0.1  s.  This  waveform  is  plotted  in  Figure  5.1.  A value 
of  y * 50  was  used  in  this  example  since  after  fifty  time  steps  the 
waveform  has  essentially  damped  to  zero.  Table  5.2  shows  the  results 
of  several  statistical  tests  on  these  data  using  a standard  deviation 
of  the  noise  ranging  from  0.04  down  to  0.001.  It  was  found  that  if  a 
is  greater  than  0.04,  a signal-to-noise  level  of  3.3  dB,  the  expected 
values  of  the  two  returned  poles  were  not  even  complex  conjugates  of 
each  other.  Thus,  when  the  signal-to-noise  level  was  worse  than  3.3  dB , 
Prony's  method  was  completely  corrupted  by  the  noise.  Note  in  Table  5.2 
that  as  the  signal-to-noise  level  gets  better  the  standard  deviation 
of  the  poles  lowers.  It  appears  that  at  around  15.0  to  20.0  dB  the  per- 
cent error  and  the  standard  deviation  of  the  poles  are  at  a tolerable 
level,  which  is, of  course,  subject  to  the  requirements  of  the  problem 
being  studied.  Notice  in  this  example,  as  in  Example  1,  that  the 
theoretical  value  of  the  N + l*’*1  eigenvalue  compares  closely  to  the  cal- 
culated mean  value.  In  all  cases  the  theoretical  value  is  slightly 
lower  than  the  calculated  mean  value. 

Example  3 

For  this  example  the  transient  response  of  Figure  4.1  in  Chapter  4 
was  used.  The  six  pole  pairs  and  their  associated  amplitudes  are 
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TABLE  5.3 

AMPLITUDE  OF  EACH  POLE  COMPONENT  FOR  EXAMPLES  3 AND  4 


Pole 

Amplitude 
Example  3 

Amplitude 
Example  4 

-0.082  ± j 0.926 

1.0 

1.0 

-0.147  ± j 2.874 

1.0 

0.5 

-0.188  ± j 4.835 

1.0 

0.25 

-0.220  ± j 6.800 

1.0 

0.125 

-0.247  ± j 8.767 

1.0 

0.0625 

-0.270  ± j 10.733 

1.0 

0.03125 

! 


N rs  go  O rs  o 

oo  <r  00  CN 

O H H fN  fN  CN 


sO  ^ in  O r*  cn 

n n n o vO  n 

on  co  co  co  n n 


o o o o o o 

I I I I I I 


O N ^ vO  CO  O 


o 

o , 

m 

CO 

vO 

CN 

CN 

m 

o 

CN 

cn 

Os 

• 

*n 

n 

sO 

ON 

rs 

O 

in 

cn 

m 

i — 

CN 

o 

00 

<r 

CO 

— H 

<4 

f-ss 

CN 

r^» 

cn 

O 

sO 

cn 

11  i 

o 

n 

<— 1 

CN 

CN 

CN 

ON 

00 

oo 

00 

r. 

r>. 

o 

o 

o 

o 

O 

o 

o 

CN 

•<r 

vO 

00 

o 

o 

1 

i 

1 

1 

1 

i 

<r  04 

in  m cn  o O o 
o o o o o o 


CO  N CN  O 

O O o o o o 


o o o o o o 


o o o o 


os  m vO  o cn  vo 

ON  lA  (*1  N ON  o 

-<r  00  r-t  <r 

O r4  H (N  CN  CM 


CO  <r  A rs  \C 

n CN  H AN  CO  CN 

cn  rs  cn  O o cn 

cn  00  oo  co  n 


CM  in  ON  <f  CN  O 
CN  o o o o o 


n sin  h o 
CN  f-l  O O 3 3 


o o o o o o 

I I I I I I 


O CN  sj  vO  CO  O 


o o o o o o 


o o o o o 


*j  o m cn  rj  o 

r»  N CO  (N  H H 

r^.  vj-  rH  m 

O H H CN  CN  CN 


in  CN  AN  CN  H 
CN  O AN  \D  O'  CN 
CN  r*»  (N  O iO  AN 
ON  CO  CO  CO  rs  N 


ON 

cn  cn  on  h in  o 
<r  o rH  h o o 


A'  AN  vf  ON  CN  o 

cn  ^ cn  o o o 


o o o o o o 

I I I I I I 


O (N  vO  CO  O 


o o o o o o 


o o o o o 


CQ 

"0 

CO 

o 

>4- 

o 

>4 

CN 

m 

ON 

m 

r. 

<r 

s4* 

O 

r4 

• 

• 

cn 

<n 

CN 

rH 

m 

*n 

•H 

m 

s4- 

*— 1 

—4 

as 

o 

P's 

sO 

O 

in 

CN 

m 

O 

■n 

r* 

cn 

CN 

| O 

**h 

CN 

CN 

CN 

Os 

00 

CO 

CO 

rs. 

o 

o 

o 

o 

o 

O 

o 

CN 

s4- 

b 

00 

o 

D 

1 

i 

1 

1 

1 

1 

CN  ON  an  N O'  CN 
O'  in  cn  o O 


O'  cn  an  vo  m cn 

>j  H o o 


o o o o o o 


o o o o o 


in  00  f>»  AN  CN  CN 
■o  cn  vo  <r  O'  a4 
cn  cn  <r  on  m r* 

H H H H CN  CN 


O ^ H O UN  (N 

m cc  cn  <r  cn  o 

m co  cn  h in  n 

O'  a*  co  As 


rs  h m co  vj  cn 
in  CO  CO  AN  H o 


CN  vO  ON  O vO  CN 
H O'  CO  M O O 


o o o o o o 

I I I I I I 


O CN  sj  vO  CO  O 


<r  o o o o o 


CN  o O © O 


u 

A3  '4-1 

c o 

yj 

30  u OJ 
fl  u H 
6 A3  O 

»-*  0-  Oi 


u 

•u  f3 

C B- 

Z)  u 

u o — < 

L4  Li  A3 
0)  Li  1) 

a.  uj  oi 


V 

U 

•H 

u 

0 

30 

4—1 

u 

u 

A3 

u 

a) 

u 

S 

A3 

CL, 

1-4 

Cl, 

SIGNAL  TO  NOISE  RATIO  OF  THE  DIFFERENT  COMPONENTS 
OF  THE  SIGNAL  IN  EXAMPLE  3 AS  A FUNCTION 
OF  THE  STANDARD  DEVIATION  OF  THE  NOISE 


TABLE  5.6 

THE  THEORETICAL  AND  EXPERIMENTAL  VALUES  OF  THE 
N 4-  1 EIGENVALUE  FOR  EXAMPLE  3 AND  EXAMPLE  4 


Example 

3 


Example 

4 


Standard 
Deviation  of 
Noise  a 

Theoretical 
Value  of 
N + 1 

Eigenvalue 

Mean  Value  of 
N + 1 

Eigenvalue 

Standard  Deviation 
of  N + 1 
Eigenvalue 

0.01 

0.02 

0.0218 

4.92  x 10-3 

0.008 

0.0128 

0.0139 

3.15  x 10"  3 

0.005 

0.0050 

0.00546 

1.23  x 10-3 

0.003 

0.0018 

0.00197 

4.43  x 10"4 

0.001 

0.2  x 10"  3 

0.22  x 10-3 

4.92  x 10-5 

0.002 

0.8  x 10_3 

0.871  x 10“3 

1.96  x 10"4 

0.001 

0.2  x 10-3 

0.218  x 10"3 

4.93  x 10‘5 

0.0005 

0.5  x in "4 

0.546  x 10"A 

1.23  x 10‘5 

BBI99H 


tabulated  in  Table  5.3.  Note  that  for  this  example  the  amplitude  of 
each  pole  is  the  same.  Noise  with  standard  deviation  from  0.01  to  0.001, 
giving  a signal-to-noise  level  from  27.4  dB  to  47.4  dB,  was  added  to  the 
signal  and  statistical  tests  were  run  on  the  data.  For  all  cases  a value 
of  y * 200  was  used  and  the  number  of  Monte  Carlo  trials  was  twenty. 

Table  5.4  shows  the  results  of  these  tests.  For  this  particular  signal, 

poles  could  not  be  extracted  for  levels  of  noise  greater  than  a = 0.01. 

That  is,  for  noise  levels  greater  than  a = 0.01,  poles  were  produced 
but  they  did  not  appear  in  complex  conjugate  pairs  and  had  no  apparent 
relationship  to  the  original  poles. 

Since  the  total  waveform  was  made  up  of  the  sum  of  six  damped 
sinusoids,  it  is  useful  to  study  the  signal-to-noise  ratio  for  each  of 
the  six  individual  components.  Table  5.5  shows  that  for  the  case  of 
a = 0.01  with  a total  signal-to-noise  ratio  of  27.4  dB  the  highest 
frequency  component  has  a signal-to-noise  ratio  of  only  14.8  dB.  This 
level  is  about  the  same  as  the  level  in  Examples  1 and  2 where  the  poles 

could  be  extracted  at  a signal-to-noise  level  of  about  2 to  3 dB,  but 

in  this  example  the  highest  signal-to-noise  level  of  an  individual  com- 
ponent is,  as  stated  above,  14.8  dB.  Hence,  indications  are  that  sum- 
ming several  signals  together  in  effect  makes  it  necessary  to  have  a 
better  signal-to-noise  ratio  to  allow  extraction  of  the  poles.  This 
point  is  further  demonstrated  in  Example  4. 

The  values  of  the  N + 1^  eigenvalues  calculated  in  Example  3 are 
tabulated  in  Table  5.6.  Note  here,  as  in  the  previous  two  examples, 
that  the  mean  value  of  the  eigenvalues  is  very  close  but  always  slightly 
higher  than  the  theoretically  predicted  values.  In  all  cases  it  was 
found  that  the  N*"*1  eigenvalue  was  sufficiently  distinct  from  the  N + lt'1 
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eigenvalue  so  that  proper  identification  of  the  cutoff  point  could  be 

made.  For  example  in  the  worst  case,  that  of  a * 0.01,  the  N + 1th 

t h 

eigenvalue  has  a mean  value  of  0.0218  and  the  N eigenvalue  has  a mean 
value  of  0.070.  Thus,  if  the  standard  deviation  of  the  noise  is 
known  in  advance,  there  should  be  no  problem  in  detecting  the  cutoff 
point  which  gives  the  value  of  N. 

Example  A 

For  the  final  example  the  same  six  pole  pairs  of  the  previous  example 
were  used  but  the  amplitudes  of  the  poles  were  reduced  as  indicated  in 
Table  5.3.  The  resulting  transient  response  curve  is  shown  in  Figure  5.2. 
Here  again  a value  of  y = 200  was  used  and  twenty  Monte  Carlo  trials  were 
performed.  Table  5.7  shows  the  results  of  noise  with  a * 0.002,  0.001, 
and  0.0005.  The  noise  level  at  which  the  true  poles  could  be  first  de- 
tected for  this  example  is  o * 0.002  which  corresponds  to  a signal-to- 
noise  level  of  39.7  dB.  This  level  is  far  higher  than  any  needed  in  the 
previous  examples,  and  is  due  to  the  fact  that  some  of  the  components 
of  the  signal  have  very  small  signal  levels.  Table  5.8  gives  the  signal- 
to-noise  level  for  each  of  the  six  components  in  this  example.  Note 
that  the  highest  frequency  has  a signal-to-noise  level  of  -1.3  dB  which 
indicates  that  the  noise  level  is  higher  than  the  average  signal  level. 
However,  the  signal  must  have  been  above  the  noise  level  for  a significant 
portion  of  the  data  window  in  order  for  it  to  be  detected.  The  standard 
deviations  of  the  extracted  poles  for  the  o = .002  noise  level  indicate 
that  the  method  is  not  accurate  enough  at  this  level  and  indeed  does  not 
appear  to  be  accurate  until  a noise  level  of  c ■ 0.0005  is  reached.  At 
this  level  the  highest  signal-to-noise  level  for  any  single  component 
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TABLE  5.7 


VALUES  OF  EXTRACTED  POLES  WITH  THEIR  PERCENT  ERROR  AND 
STANDARD  DEVIATION  AS  A FUNCTION  OF 
SIGNAL  TO  NOISE  RATIO  FOR  EXAMPLE  4 


Signal  Co  Noise  Racio 


Standard  Deviation 
of  Noise  a 


Real  Part  of  Poles 


Imaginary  Part  of  Poles 


Percent  Error 
Real  Part 


Percent  Error 
Imaginary  Part 


Standard  Deviation  of 
Real  Part  of  Poles 


Standard  Deviation  of 
Imaginary  Part  of 
Poles 


39.7  dB 


-0.0809 

-0.1493 

-0.3629 

-0.1483 

-0.3109 

-0.2315 


0.9220 

2.8675 

4.9708 

6.7206 

8.6936 

10.6593 


0.0113 

0.1034 

0.8509 

0.1000 

0.4254 

0.0834 


0.0310 

0.0382 

0.7797 

0.4604 

0.3860 

0.2846 


45.7  dB 


-0.0814 

-0.1452 

-0.1793 

-0.2036 

-0.2823 

-0.2700 


0.9251 

2.8726 

4.8293 

6.8173 

8.7688 

10.7108 


0.0054 

0.0415 

0.0600 

0.0371 

0.1871 

0.0470 


0.0120 

0.0120 

0.0585 

0.1237 

0.0632 

0.1225 


51.7  dB 


0.0005 


-0.0817 

-0.1459 

-0.1848 

-0.2159 

-0.2565 

-0.2719 


0.9256 

2.8734 

4.8342 

6.8065 

8.7686 

10.7271 


0.0026 

0.0200 

0.0300 

0.0175 

0.0785 

0.0231 


0.0057 

0.0055 

0.0283 

0.0616 

0.0211 

0.0535 


TABLE  5.8 


SIGNAL  TO  NOISE  RATIO  OF  THE  DIFFERENT  COMPONENTS 
OF  THE  SIGNAL  IN  EXAMPLE  4 AS  A FUNCTION 
OF  THE  STANDARD  DEVIATION  OF  THE  NOISE 


Portion  of  Signal 
Due  to  Particular 
Pole 


Standard 
Deviation  of 
Noise  a 


Sum  of  Below 
-0.082  ± j 0.926 
-0.147  ± i 2.874 
-0.188  t j 4.835 
-0.220  ± j 6.800 
-0.247  ± j 8.767 
-0.270  ± j 10.733 


0.002 

0.001 

0.0005 

Signal  to 

Noise  Ratio  dB 

39.7 

45.7 

51.7 

38.9 

44.9 

50.9 

28.1 

34.1 

40.1 

19.9 

25.9 

31.9 

17.9 

23.9 

29.9 

5.5 

11.5 

17.5 

-1.3 

4.7 

10.7 

is  better  than  10  dB.  Table  5.6  shows  that  here  again  the  mean  value 
of  the  lowest  eigenvalue  is  always  slightly  higher  than  the  theoretical 
value. 

Summary  of  Numerical  Examples 

Four  different  numerical  examples  were  studied  to  determine  the 

effects  which  noise  has  on  Prony's  method  and  to  see  if  the  value  of  N 

could  be  found  from  the  eigenvalues.  In  all  cases  the  mean  value  of 

s t 

the  calculated  N + 1 eigenvalue  was  slightly  higher  than  the  theoretical 
2 

value  of  y a . It  was  also  found  that  the  N + 1 eigenvalue  was  suf- 
ficiently different  from  all  others  of  lower  order  so  that  it  is  possible 
to  identify  it.  For  all  four  examples  a different  highest  value  of 
tolerable  signal-to-noise  level  was  obtained.  However,  in  all  the 
examples  the  best  results  were  obtained  when  the  signal-to-noise  level 
was  between  10  and  20  dB  for  the  lowest  level  component  of  the  signal. 

The  examples  Indicate  that  different  tolerable  noise  levels  will 
be  found  for  different  shaped  waveforms.  One  definite  conclusion  is 
that  if  the  noise  level  is  higher  than  the  signal  level  over  the  entire 
data  window  then  that  signal  or  its  accompanying  pole  will  not  be  de- 
tectable. These  examples  also  indicate  that  further  study  should  be 
done  with  experimentally  obtained  signals  containing  truly  random  noise. 


6.  APPLICATIONS 


Now  that  the  techniques  for  properly  applying  Prony's  method  have 
been  developed  it  is  of  interest  to  study  some  of  the  electromagnetic 
applications  of  this  method.  There  are  four  major  areas  to  which  the 
extraction  of  poles  from  a transient  electromagnetic  signal  can  be  of 
great  benefit.  These  are:  system  analysis,  radar  target  recognition, 

the  study  of  spectral  characteristics,  and  data  reduction  and  extrapolation. 
These  four  applications  are  discussed  in  some  detail  in  this  chapter. 

6.1  System  Analysis 

6.1.1  Response  to  various  exciting  waveforms 

In  circuit  theory  after  the  impulse  response  is  known,  the  response 
of  the  circuit  to  any  given  driving  function  can  be  determined  by  mul- 
tiplying the  Laplace  transform  of  the  driving  function  by  the  sum  of 
poles  and  their  corresponding  residues  of  the  impulse  response.  That  is, 

N A 

R(s)  = F(s)  V i—  (6.1) 

i-1  " i 

where  R(s)  is  the  resulting  response  function  and  F(s)  is  the  Laplace 
transform  of  the  driving  function.  In  Chapter  2 it  was  shown  that  a 
similar  relation  is  true  for  electromagnetic  structures.  The  differences 
are  that  since  antennas  and  scatterers  are  distributed  systems  the  im- 
pulse response  is  a function  of  position  on  the  structure  and  a function 
of  the  incident  angle  of  the  driving  function.  Note  that  the  pole  portion 
of  the  impulse  response  is  independent  of  position  but  that  the  natural 
modes  are  functions  of  position  and  the  coupling  coefficients  are  functions 
of  the  incident  angle.  Thus,  the  general  response  function  for  an  antenna 
or  scatterer  can  be  expressed  as 
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R(s,r,p)  = F(s,r,p) I £ - _ ■- — I (6.2) 

where  the  coupling  coefficients  n^(s^,p)  and  the  natural  modes  v^(r) 
have  been  combined  into  the  one  term,  A^r.p).  Since  the  poles  s^  are 
invariant  of  position,  they  may  be  determined  by  studying  the  impulse 
response  at  any  position  on  the  structure.  After  the  poles  have  been 
evaluated,  the  residues  A^  must  be  determined  at  every  desired  point 
on  the  structure.  This  however  is  not  a difficult  problem.  Thus,  the 
poles  of  a structure  need  only  be  extracted  once  and  the  residues  must 
be  determined  at  each  position  that  the  response  is  desired  and  for 
each  angle  of  incidence  that  the  driving  function  will  be  applied.  As 
an  illustration  consider  the  following  example. 

Consider  that  the  induced  current  at  several  positions,  e.g.,  M, 
on  a dipole  due  to  several  different  time-varying  broadside  incident 
plane  waves  is  desired.  The  first  step  is  to  determine  the  induced 
current  at  all  M desired  positions  on  the  structure  resulting  from  a 
broadside  incident  impulse  plane  wave.  It  will  actually  be  necessary 
to  use  a narrow  Gaussian  plane  wave  as  an  approximation  to  the  impulse. 

The  true  impulse  response  can  then  be  obtained  by  using  one  of  the 
methods  outlined  in  Section  2.2.  After  the  impulse  response  of  the 
induced  current  is  obtained,  then  the  poles  of  the  structure  can  be  ob- 
tained by  applying  Prony's  method  to  the  current  at  one  of  the  h points 
on  the  structure.  Now  that  the  poles  have  been  determined,  the  M 
sets  of  residues  are  calculated  for  the  M positions  on  the  antenna. 
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thus  giving  one  set  of  N poles  and  M sets  of  N residues  A^r.p).  The 
induced  current  for  an  arbitrary  broadside  incident  waveform  is  then 
obtained  by  applying  (6.2).  If  an  incident  angle  other  than  broadside 
is  required,  then  it  is  necessary  to  recalculate  the  M sets  of  residues 
for  an  incident  impulse  plane  wave  at  the  new  incidence  angle.  This 
seems  a horrendous  amount  of  calculations  but  when  compared  with  the 
amount  required  to  completely  resolve  the  problem  for  each  incident  wave 
shape  it  is  really  a very  small  amount. 

6.1.2  Compatibility  with  circuit  theory 

In  practice  antennas  are  almost  always  coupled  to  some  sort  of 
circuit  or  network.  This  presents  a problem  to  the  circuit  designer 
because  the  antenna  port  is  not  classified  as  a lumped  element.  This 
problem  is  usually  not  severe  since  the  network  and  antenna  are  gen- 
erally designed  to  operate  at  one  given  frequency.  In  this  case  the 
antenna  port  can  be  characterized  by  a lumped  impedance  at  that  one 
frequency.  If,  however,  the  antenna  is  to  be  operated  over  a broad  band 
of  frequencies  the  circuit  designer  must  know  the  antenna's  terminal 
characteristics  over  that  entire  band.  The  characteristics  are  normally 
obtained  from  a set  of  graphs  relating  the  real  and  imaginary  parts  of 
the  input  admittance  to  the  frequency.  Obviously,  this  is  not  an  ideal 
situation.  The  circuit  designer  would  like  to  have  the  antenna's  input 
impedance  given  to  him  in  the  Laplace  transform  domain.  This  is  a pos- 
sibility if  Prony's  method  is  used.  All  that  is  necessary  is  to  apply 
a time-varying  voltage  to  the  antenna  terminals  and  to  measure  the  re- 
sulting current  as  a function  of  time.  If  Prony's  method  is  applied  to 
the  current  and  a set  of  poles  and  residues  is  obtained,  then  the  input 
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admittance  can  be  expressed  as 
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where  V(s)  is  the  Laplace  transform  of  the  applied  voltage.  Expression 
(6.3)  is  obviously  much  easier  for  a circuit  designer  to  use  than  a 
set  of  complicated  graphs  or  tables. 


6.1.3  Study  of  system  parameters 

It  is  quite  useful  to  know  what  effect  the  loading  of  a structure 
has  on  the  positions  of  its  poles  in  the  complex  plane.  Likewise  it 
is  beneficial  to  know  what  effect  the  positions  of  the  poles  have  on 
the  behavior  of  the  transient  response.  If  trajectories  of  the  poles 
could  be  determined  from  a few  experimental  cases,  it  would  be  possible 
to  predict  pole  locations  as  a function  of  loading.  Once  the  pole  lo- 
cations are  known  then  the  transient  waveforms  can  be  constructed.  One 
practical  example  of  this  problem  is  the  resistive  loading  of  a linear 
antenna  in  order  to  produce  a transmitted  field  which  simulates  an 
electromagnetic  pulse.  Tesche  [23]  has  approached  this  problem  by  de- 
termining the  pole  locations  of  a dipole  which  is  loaded  with  uniform 
resistive  loading  along  the  antenna.  He  produced  the  trajectories  of 
the  poles  by  the  "classical"  frequency  domain  search  procedure.  An 
alternative  approach  is  to  obtain  the  time-domain  solutions  for  the 
induced  currents  on  and  the  scattered  fields  from  several  uniformly 
resistively  loaded  dipoles.  Prony's  method  is  applied  to  either  the 
transient  response  of  the  induced  currents  or  the  scattered  fields  in 
order  to  determine  the  location  of  the  poles.  The  results  of  such  a 
procedure  are  given  in  the  following  example. 
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A 1.0  m dipole  with  a half-length-to-radius  ratio  of  100  was 
numerically  modeled  and  excited  by  a broadside  incident  Gaussian  pulse. 

The  induced  current  at  the  center  of  the  antenna  and  the  back  scattered 

♦ 

electric  field  was  calculated  as  a function  of  time  for  uniform  resistive 
loading  of  0,  125,  250^,  500,  750  and  1683  0/m.  These  currents  and  fields 
are  shown  in  Figures  6.1  and  6.2,  respectively.  The  time  step  size  used 
was  At  = 6.9444  x 10  ^ s.  The  loading  value  of  1683  0/m  was  chosen 
because  Tesche  [23]  calculated  that  1683  0/m  is  the  value  at  which  the 
dipole  would  become  critically  damped,  giving  a double  pole  on  the  neg- 
ative real  axis.  The  poles  that  were  extracted  using  Prony's  method  are 
plotted  in  Figure  6.3.  The  trajectories  of  the  first  seven  even  poles 
are  shown.  Note  that  for  the  value  of  1683  0/m  the  first  pole  has  split 
and  moved  toward  the  origin  and  toward  infinity,  which  indicates  that  this 
value  of  loading  does  not  give  a critically  damped  dipole  but  produces 
an  overdamped  situation.  This  does  not  imply  that  Tesche1 s calculated 
value  of  1683  0/m  is  incorrect.  All  that  is  indicated  is  that  this  nu- 
merical modeling  procedure  and  Tesche 's  numerical  modeling  procedure 
differ.  It  should  be  noted,  however,  that  the  point  at  which  the  pre- 
dicted trajectory  of  the  first  pole  in  Figure  6.3  splits  is  approximately 
the  same  value  which  Tesche  indicates  in  his  paper. 

One  very  interesting  point  in  the  above  example  is  that  for  even 
the  overdamped  transient  response  Prony's  method  was  capable  of  ex- 
tracting the  first  five  poles  of  the  system.  This  is  a very  important 
point  since  if  the  1683  fl/m  case  in  Figure  6.1  is  studied  it  appears  that  no 
oscillations  occur.  Yet  Prony's  method  was  capable  of  determining  the 
first  five  modes  including  the  splitting  of  the  first  pole.  This  example 
then  indicates  that  Prony's  method  will  work  with  heavily  loaded  structures 
and  with  fat  or  thick  structures. 
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Figure  6.1.  Induced  current  at  center  element  of  1 m dipole 
with  resistive  loadings. 
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6. 2 Radar  Target  Recognition 

The  problem  of  electromagnetic  recognition  of  a radar  target  has 
received  a great  deal  of  attention  in  recent  years.  To  recognize  a 
target  with  a single-frequency  radar  requires,  in  principle,  bistatic 
scattering  cross-section  information  for  a single  transmitter  aspect 
angle  and  for  4tt  steradians.  In  practice,  with  a priori  information 
about  the  body,  it  is  feasible  to  match  a body  to  one  of  a known  group 
of  bodies  using  a somewhat  smaller  angular  range.  With  multiple- frequency 
monostatic  radar  cross-section  information  available  over  a limited 
angular  range,  it  is  possible  to  identify  a body  as  well  [29].  Theoret- 
ically, the  monostatic  reconstruction  requires  information  over  Att 
steradians  and  an  infinite  frequency  range. 

Recent  work  indicates  the  feasibility  of  using  short-pulse  (impulse 
response)  excitation  in  a monostatic  configuration  for  purposes  of  target 
recognition  from  a class  of  bodies  [30]  - [31]. 

In  the  first  of  these  schemes,  reported  by  Sperry  Research  Center, 
information  about  the  body  is  contained  in  the  time-history  of  the 
backscattered  waveform,  i.e.,  a time  signature.  This  technique  is  com- 
plicated by  the  fact  that  the  time  signature  is  dependent  on  the  aspect 
angle  of  the  target.  The  result  is  the  requirement  of  a rather  large 
catalog  of  known  signatures. 

Mains  and  Moffatt  [31]  introduced  the  concept  of  using  the  complex 
natural  resonances  of  a target  as  a basis  for  target  recognition.  They 
make  use  of  the  fact  that  a few  natural  resonances  of  a body  are  adequate 
to  distinguish  the  body  within  a finite  collection  of  bodies.  They  also 
make  use  of  the  fact  that  the  natural  resonances  of  a body  as  manifested 
in  a scattered  waveform  are  aspect  independent:  they  do  not  depend  on 


116 


the  angular  orientation  of  the  target.  They  justify  this  fact  empirically 
in  [31].  Baum,  who  introduced  the  Singularity  Expansion  Method  (SEM) 

[1]  - [2],  has  rigorously  demonstrated  this  aspect  independence  in  the 
context  of  the  SEM  method. 

Thus,  this  method  provides  a convenient  means  of  1)  characterizing 
a transient  scattered  waveform  in  terms  of  a few  coefficients  and  complex 
resonant  frequencies  in  a series  of  decaying  sinusoids,  and  2)  separating 
aspect-dependent  characteristics  of  the  waveform  from  characteristics 
intrinsic  to  the  body.  In  particular,  the  complex  frequencies  of  the 
sinusoids  in  the  series  are  intrinsic  to  the  body  while  the  amplitude 
of  the  sinusoids  depends  on  the  wave  shape  of  the  excitation  and  on  the 
aspect  angle. 

Two  significant  practical  features  accrue  from  the  use  of  the  SEM 
representation  with  primary  attention  directed  to  the  complex  frequency 
qualities:  1)  aspect  dependence  is  suppressed  and  2)  the  incident 

wave  shape  need  not  be  impulse  - it  needs  only  to  have  sufficient  spectral 
content  to  excite  several  of  the  complex  sinusoidal  modes. 

Since  Prony's  method  can  be  used  to  numerically  extract  the  weighting 
coefficients  and  the  complex  frequencies  from  a digitized  backscattered 
time-signature,  Pearson,  Van  Blaricum  and  Mittra  [32]  suggested  that 
the  method  be  used  to  aid  in  the  target  recognition  problem.  The  ex- 
tracted frequencies  (rather  than  the  complete  time-signature)  would 
serve  as  the  input  to  a pattern-recognition  algorithm.  This  method  is 
believed  to  represent  a significant  improvement  over  direct  time-history 
pattern  recognition  because  the  aspect-dependence  is  suppressed  and 
because  the  frequencies  comprise  a smaller  set  of  numbers  than  the 
entire  sampled  time-signature. 
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A block  diagram  for  the  recognition  procedure  is  given  in  Figure 
6.4.  The  system  consists  of  a transmitter,  a receiver  capable  of 
digitizing  the  transient  return,  and  a digital  processor  with  some  con- 
venient operator  communications  device  such  as  a CRT.  Of  course,  a 
tracking  system  and  an  antenna/duplexer  are  present. 

The  waveform  generator  must  produce  a transient  waveform  with 
reasonably  broad  spectral  content.  The  power  and  frequency  requirements 
are,  of  course,  dictated  by  the  application.  In  particular,  target 
range  requirements  and  receiver  performance  determine  a power  specifi- 
cation, while  the  body  size  determines  frequency.  The  frequency  must 
be  such  that  f = c/D  with  D a characteristic  dimension  of  the  class  of 
targets  and  c the  velocity  of  light.  An  alternative  to  a transient 
source  and  receiver  is  a multiple-frequency  pulsed  carrier  scheme  such 
as  Mains  and  Moffatt  suggest  [31].  This  scheme  uses  CW  scattering  re- 
turns to  synthesize  the  transient  response  of  the  body. 

The  transient  digitizer  provides  a means  of  detecting,  sampling, 
and  digitally  expressing  the  transient  return  signal.  In  its  simplest 
form,  this  element  might  be  a storage  oscilloscope  with  A-D  converters 
on  its  vertical  and  horizontal  drive  signals.  There  are  other  more 
sophisticated  systems  available  commercially. 

The  digital  processor  must  be  able  to  perform  the  following 
operations:  a)  extract  the  poles  of  the  waveform  by  using  Prony's 

method;  b)  identify  (based  on  data  designed  into  the  system)  poles  due 
to  the  transmitted  wave  shape  and  to  the  RF  system;  and  c)  conduct  a 
pattern  recognition  procedure  to  provide  a set  of  "likelihood  votes" 
associating  the  target  with  a catalog  of  expected  targets.  The  first 
item  above  is  self-explanatory. 
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Figure  6.4.  Schematic  representation  of  the  radar  target  recognition 
svstem. 
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Any  waveform  that  can  be  radiated  is  oscillatory  in  character 
and  will  accordingly  have  associated  with  it  at  least  one  "waveform 
pole."  For  example,  a sinusoidal  segment  will  possess  a pole  at  s * ju^, 
where  Uq  is  the  radian  frequency  of  the  sinusoidal  signal.  In  addition, 
it  is  conceivable  though  undesirable  that  the  receiving  antenna  and  RF 
hardware  might  introduce  poles  into  the  waveform  representation.  All 
of  these  "system"  poles  are  known  during  system  design,  however,  and 
it  is  a straightforward  algorithm  to  locate  and  delete  them  from  the 
set  of  poles  extracted  from  the  signal.  An  alternative  approach  is  to 
numerically  deconvolve  the  waveform  spectrum  and  system  transfer  function 
from  the  received  waveform. 

The  final  process  is  to  compare,  by  means  of  some  pattern  analysis 
algorithm,  the  observed  pole  set  with  a data  base  of  known  targets' 
poles.  This  algorithm  would  display,  as  its  output,  a "vote"  or  prob- 
ability for  a known  target  or  targets.  This  information  could  be  sup- 
plemented by  an  indication  of  recognition  confidence,  for  example,  the 
number  of  poles  successfully  used  in  the  identification.  A second  sup- 
plementary output  might  be  tracking  data  on  the  target. 

For  this  method  to  be  successful  as  a target  identification  scheme, 
a few  poles  (say  less  than  ten)  must  characterize  a given  target  among 
all  potential  targets.  Mains  and  Moffatt  [31]  discuss  this  problem  in 
some  detail  and  present  the  pole  configuration  for  several  thin-wire 
geometries  to  indicate  that  a few  poles  do  distinguish  among  similar 
objects.  Also,  another  criterion  for  this  system  to  be  successful  is 
a satisfactory  signal-to-noise  level  of  the  returned  time-signature  so 
that  Prony's  method  can  properly  extract  the  poles.  Bern!  [33]  has 
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suggested  an  identification  technique,  which  is  also  based  on  Prony's 
algorithm,  although  Prony's  method,  as  such,  is  not  mentioned  in  his 
paper.  His  approach  is  to  use  a series  of  exciting  pulses  which  then 
produces  a series  of  response  functions  which  are  correlated  to  produce 

I 

what  he  calls  a "target  correlation  matrix."  His  correlation  matrix 
is  similar  to  the  $ matrix  of  Prony's  method  of  Section  4.2.  The  dif- 
ference is  that  his  correlation  matrix  is  filled  with  samples  from 
several  different  sets  of  transient  response  data  while  that  of  Prony's 
method  uses  just  one  set  of  transient  response  data.  Bemi's  approach 
may  thus  prove  useful  in  reducing  the  noise  level  and  some  of  the  signal- 
to-noise  ratio  problems  outlined  in  Section  6.3.  It  appears  that  the 
use  of  Prony's  method  as  a tool  for  radar  target  recognition  may  be  of 
great  importance  but  more  study  of  the  problem  is  needed  in  order  to 
satisfy  many  of  the  unanswered  questions. 

6. 3 Study  of  Spectral  Characteristics 

It  has  been  pointed  out  already  that  after  the  poles  and  residues 
of  a system  have  been  determined  it  is  then  possible  to  write  the  fre- 
quency domain  version  of  the  system's  response  as 
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where  the  poles  s^  have  been  written  in  terms  of  their  real  and  imaginary 
parts  as 

si  m + 1 (6.5) 


Thus,  the  frequency  domain  response  can  be  obtained  directly  from  the 
time  domain  response  without  having  to  perform  a Fourier  transform. 
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The  real  advantage  is  that  in  order  to  perform  a Fourier  transform 
accurately  it  is  necessary  to  have  a very  long  time  history  of  the 
transient  response.  That  is,  it  is  usually  necessary  to  have  enough 
time  samples  so  that  the  transient  response  has  decayed  to  zero  or  to 
its  steady-state  value.  The  examples  of  Section  3.3  showed  that  Prony's 
method  is  capable  of  extracting  the  poles  from  a time  window  that  is 
extremely  short  when  compared  to  the  whole  response..  Thus,  if  it  is 
too  expensive  to  calculate  or  measure  more  than  just  a few  transient 
data  samples, it  is  generally  impossible  to  perform  an  accurate  Fourier 
transform,  but  if  Prony's  algorithm  is  used,  it  is  then  possible  to 
extract  the  spectral  characteristics  from  this  narrow  time  window. 

If  the  transient  response  is  not  the  impulse  response  and  the  im- 
pulse response  or  the  frequency  domain  transfer  function  is  desired, 
it  is  possible  to  perform  the  deconvolution  after  the  poles  have  been 
found.  For  example,  if  Prony's  method  gives  a frequency  domain  response 
function  of  the  form  of  (6.4)  and  if  the  frequency  domain  response  of 
the  driving  function  F(s)  is  known,  it  is  possible  to  obtain  the  fre- 
quency domain  transfer  function: 


H(s) 


F (s) 


(6.5) 


Note  the  similarities  between  (6.5)  and  (6.3),  the  expression  for  the 
input  admittance  of  an  antenna.  If  F(s)  is  the  applied  voltage  at  the 
driving  point  and  the  and  the  s^  are  the  residues  and  poles  for  the 
induced  current  at  the  driving  point,  then  H(s)  is  just  the  driving- 
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parts,  respectively,  of  the  input  admittance  of  the  1.0  m dipole  discussed 
in  Section  3.3.  The  solid  line  represents  the  admittance  obtained  by 
using  the  conventional  Fast  Fourier  Transform,  and  the  dotted  line  is 
the  admittance  obtained  using  Expression  (6.3).  Since  the  driving 
function  was  a Gaussian  pulse,  it  was  possible  to  express  its  Laplace 
transform  F(s)  analytically.  Note  that  the  two  methods  give  results 
that  compare  closely  for  all  but  the  higher  frequencies. 

Similarly,  if  the  normalized  radar  cross  section  of  a scatterer 
is  desired,  it  can  be  obtained  by  using  the  expression 
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where  the  dependence  on  the  incidence  and  observation  angles  and  polar- 
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ization  has  been  suppressed.  The  E (w)  term  is  simply  the  frequency 
domain  response  of  the  applied  incident  field  and  Erad(“)  is  the  frequency 
domain  response  of  the  scattered  field  measured  at  distance  r from  the 
scatterer.  E ^ (id)  can  be  obtained  using  Prony's  method  on  the  calculated 
or  measured  scattered  transient  response  data. 

When  a conventional  Fourier  transform  routine  is  used,  the  result 
is  always  a set  of  frequency  domain  data  points  at  a specified  frequency 
interval  which  is  dependent  on  the  time  step  size  of  the  transient  data 
used.  Many  times  the  study  of  the  spectral  characteristics  at  just  a 
few  select  frequencies  is  of  interest.  If  any  of  the  analytical 
Expressions  (6.3),  (6.4),  (6.5),  or  (6.6)  are  used.it  is  possible  to 
calculate  the  spectral  characteristics  at  any  frequency  desired. 

Thus,  if  Prony's  method  is  used,  it  is  possible  to  obtain  as  few  or  as 
many  frequency  domain  samples  at  any  frequency  of  interest,  within  the 
bandwidth  of  the  system. 
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From  the  discussions  and  examples  mentioned  earlier,  it  is  obvious 
that  Prony's  method  takes  a set  of  transient  response  data  and  reduces 
it  to  a small  set  of  poles  and  residues.  These  poles  and  residues  can 


then  be  used  to  reconstruct  the  transient  response  for  all  time,  in- 
cluding those  time  steps  past  the  end  of  the  data  window  used.  Hence, 
Prony's  method  is  a valuable  method  for  data  reduction  and  extrapolation. 
Anyone  who  has  ever  generated  transient  response  waveforms  by  numerical 
or  experimental  procedures  soon  learns  that  the  storage  problem  required 
in  keeping  all  of  the  data  is  very  large.  If  Prony's  method  is  used, 
it  is  not  only  possible  to  reduce  the  amount  of  data  storage  required 
but  it  is  also  possible  to  reduce  the  amount  of  transient  data  produced 
in  the  first  place.  The  following  example  will  demonstrate  this  adequately. 

Consider  again  the  experimentally  produced  data  which  were  discussed 
in  Section  3.3,  Example  9.  The  experimental  transient  range  produced 
512  data  samples  over  the  time  window  shown  in  Figure  3.18.  Prony's 
method  was  then  applied  to  only  100  of  these  time  samples  at  every  fifth 
time  step  to  produce  the  forty-one  poles  shown  in  Figure  3.19.  The  forty- 
one  poles  and  residues  were  then  used  to  reproduce  the  measured  tran- 
sient response  and  to  extrapolate  the  measured  response  to  very  late 
time,  100  ns.  This  reproduced  response  is  shown  in  Figure  6.6.  Thus, 
using  only  eighty-two  complex  numbers,  the  original  measured  response 
was  reproduced  and  extrapolated  to  add  additional  information.  Remember 
also  that  only  100  of  the  original  measured  samples  were  needed. 

Hence,  a response  can  be  measured  for  a fairly  short  period  of  time; 
then, Prony's  method  can  be  applied  to  the  response  to  reduce  it  to  a 
set  of  poles  and  residues.  The  poles  and  residues  can  then  be  used  to 
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reconstruct  the  original  data,  extrapolate  the  original  data  to  later 
time  values,  and  produce  the  spectral  content  by  using  the  methods  of 
Section  6.3. 

6. 5 Other  Applications  of  Prony's  Method 

This  chapter  dealt  with  applications  of  Prony's  method  for  extracting 
the  singularities  of  a system  from  its  transient  response.  There  are 
many  other  uses  for  Prony's  method  which  have  been  developed  in  the  past. 
These  uses  range  from  representation  of  electrocardiograms  [17]  to  the 
measuring  of  the  vertical  angles  of  arrival  of  HF  sky-wave  signals  [34]. 
These  applications  are  very  interesting  and  for  that  reason  an  additional 
bibliography  is  provided  which  lists  many  of  the  papers  in  which  Prony's 
method  has  been  used.  The  techniques  developed  in  this  thesis  could  be 
used  to  eliminate  many  of  the  problems  discussed  in  the  papers. 


7.  AN  ALTERNATIVE  TO  PRONY'S  METHOD 
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In  previous  chapters,  Prony's  method  has  been  the  only  technique 
discussed  for  extracting  the  poles  from  a set  of  transient  response 
data.  However,  another  approach,  known  as  the  Pade  approximation  method 
[35],  can  be  used.  The  Pade  method  will  be  discussed  in  this  chapter 
and  will  be  shown  to  be  extremely  limited  in  its  usefulness. 

7.1  Padg  Approximation  Method 

If  the  Laplace  transform  R(s)  of  the  transient  response  function 
R(t)  is  analytic  at  the  origin,  then  R(s)  can  be  represented  in  a Taylor's 
series  expansion  about  the  origin.  The  Taylor's  series  can  be  denoted 
by 

CO 

R(s)  - rQ  + r^s  + r2s“  + . . . + rksk  + £ r^1  . (7.1) 

i=k+l 


A rational  function  in  s can  always  be  found  such  that  its  Taylor's 
expansion  has  the  same  leading  terms  as  those  of  (7.1).  This  rational 
function  is  known  as  the  Pade  approximant  of  R(s).  The  Pad£  approximant 
is  usually  written  in  the  form 
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Note  that  if  the  transient  response  contains  exactly  N poles  then 
Expression  (7.2)  is  an  analytical  expression  for  the  Laplace  transform 
of  the  sum  of  exponentials 
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Note  also  that  if  the  coefficients  of  the  denominator  of  (7.2)  can  be 
obtained  then  it  is  a simple  procedure  to  find  the  poles  of  the  system. 


All  that  would  need  to  be  done  is  to  find  the  roots  of  the  NCn  order 
polynomial 


bg  + b^s  + b^s  + 


V “i-1  (s  “ si>  = 0 • 


The  coefficients  b^  of  (7.3)  can  be  obtained  in  the  Pad£  approximant 
sense  by  setting  the  first  k + 1 terms  of  (7.1)  equal  to  (7.2).  That 
is,  let 

, . 2 , , m 

- , a.  + a,  s + a»s  + . . . + a s 

r0  + ris  + r2s  + . . . + rks  = ~ 2 N (7.5) 

bg  + b^s  + b^s  + . . . + b^s 

where  the  value  of  k must  be  equal  to  m + N.  If  the  denominator  of  the 
right  side  of  (7.5)  is  multiplied  by  the  left  term  and  like  powers  of 
s are  equated,  the  following  set  of  equations  is  obtained. 
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If  in  (7.6b)  the  bg  terms  are  set  equal  to  one,  the  N remaining  b.^  may 
be  solved  in  terms  of  the  r^.  Then  the  a^  in  (7.6a)  can  be  solved  in 
terms  of  the  b^  and  the  r^. 

The  remaining  problem  is  the  determination  of  the  r^  which  are  the 
coefficients  of  the  Taylor's  expansion  of  the  transient  response.  These 
can  be  obtained  if  the  transient  waveform  is  modeled  with  a set  of  step 
functions  u(t  - t^,)  and  a set  of  ramp  functions,  u ^(t  - t^).  The  step 
function  and  the  ramp  function  have  Laplace  transforms  which  are  known. 
These  are 

. -st . 

L[u(t  - t.)]  = ± e 1 (7.7a) 

1 -sti 

l[u_1(t  - t±)]  = -2  e (7.7b) 
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where  the  time  t^  is  the  turn-on  time  for  the  step  and  the  ramp  functions. 

Hence,  if  the  transient  response  can  be  modeled  with  a set  of  step  and 

ramp  functions,  then  its  Laplace  transform  can  be  expressed  as  a sum  of 
the  terms  (7.7a)  and  (7.7b).  As  a simple  example  of  the  modeling  pro- 
cedure consider  the  portion  of  a transient  response  shown  as  a dotted 
line  in  Figure  7.1.  This  response  has  been  modeled  as  a sum  of  step 
and  ramp  functions  as  indicated  by  the  solid  lines.  The  analytical  time 
domain  expression  for  the  model  can  then  be  written  for  this  example  as 

R(t)  - 13  u_x(t)  - 13  u_x(t  - 1)  + 13  u (t  - 1) 

+ 3 u_1(t  - 1)  - 3 u_1 (t  - 2)  - 8 u_1(t  - 2)  + . . . . (7.8) 


Expression  (7.8)  then  has  a Laplace  transform  of 
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which  reduces  to 
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The  exponential  terms  in  (7.9)  which  result  from  the  time  shift  have  a 
Taylor  series  expansion  of 


-t  s (-t  s)2  (-t  s)3 

e = 1 + (-t . s)  + tt + :r; + . . . + 


(7.10) 


If  this  expansion  is  substituted  into  the  expression  for  the  Laplace 
transform  of  the  pulse  and  ramp  models  of  the  transient  response,  then 
the  response  can  be  written  in  the  form  of  (7,1).  This  expansion  should 
then  be  truncated  to  include  only  k + 1 terms  as  needed  to  solve  the 
set  of  Equations  (7.6b).  It  is  now  a simple  procedure  to  obtain  the 
unknown  coefficients  b^  by  solving  (7.6b).  The  poles  are  then  obtained 
as  the  roots  of  the  polynomial  (7.4)  which  has  the  coefficients  b^. 

After  the  a^  have  been  solved  for  from  (7.6a),  residues  of  the  poles 
can  be  obtained  by  performing  a partial  fraction  expansion  on  the  rational 
function  (7.2). 

The  above  approach  for  finding  the  poles  appears  to  be  a simple 
procedure  but  it  will  be  shown  in  the  next  section  that  because  of  the 
approximations  made  this  technique  is  only  useful  for  systems  contain- 
ing a couple  of  poles. 
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7 . 2 Examples  of  the  Pade  Procedure 

Example  1 

As  the  first  example  of  the  Pad£  procedure  consider  the  response 
function 

R(t)  - e-2t  sin  irt  (7.11) 

which  has  been  discussed  previously  as  Example  2 of  Section  5.3  and 
plotted  as  Figure  5.1.  A total  of  200  samples  were  generated  and  the 
response  was  modeled  for  20,  40,  80,  100,  150,  and  200  of  these  samples 
using  the  pulse  and  ramp  approximations  discussed  in  Section  7.1.  The 
results  of  applying  the  Pad£  approximation  to  determine  the  poles  are 
presented  in  Table  7.1.  Note  that  the  results  are  not  satisfactory 
until  eighty  samples  were  used.  Figure  5.1  shows  that  the  response 
does  not  damp  to  nearly  zero  until  after  time  step  forty.  In  order  to 
get  valid  results  all  of  the  transient  response  must  be  used  until  the 
time  at  which  the  response  is  nearly  equal  to  zero.  The  results  for 
eighty  through  200  samples  are  promising  however. 

Example  2 

Here  a transient  response  was  generated  using  two  sets  of  complex 
pole  pairs.  The  response  used  is 

r(t)  - e 2t  sin  irt  + e 3t  sin  1.5  fft  . (7.12) 

The  results  of  the  Pad£  approximation  procedure  for  this  transient 
response  are  tabulated  in  Table  7.2.  Note  that  in  this  case  satisfactory 
results  were  not  obtained  until  200  samples  were  used  even  though  the 
response  has  damped  to  zero  at  around  the  fiftieth  time  step.  This  in- 
dicates that  a larger  number  of  samples  is  needed  to  resolve  the  presence 

133 


TABLE  7.1 


THE  RESULTING  POLES  USING  THE  PADE  APPROXIMATION  ON  THE  TRANSIENT 
RESPONSE  OF  EXAMPLE  1:  R(r)  =«  e-2t  sin  irt 


Number  of 

Poles  Extracted 

Percent 

Error 

Samples  Used 

Real  Part 

Imaginary  Part 

20 

-1.5009  ± j 2.1523 

13.40 

26.56 

40 

-1.9076  ± j 2.8563 

2.48 

7.66 

80 

-1.9998  ± j 3.1405 

0.004 

0.004 

100 

-2.0001  ± j 3.1414 

0.002 

0.004 

150 

-2.0001  ± j 3.1414 

0..002 

0.003 

200 

-2.0001  ± j 3.1414 

0.002 

0.003 

TABLE  7.2 

THE  RESULTING  POLES  USING  THE  PADE  APPROXIMATION  ON  THE  TRANSIENT 
RESPONSE  OF  EXAMPLE  2:  R(t)  - 2-2t  sin  irt  + e~3t  sin  1.5  irt 


Number  of 

Percent 

Error 

Samples  Used 

Poles  Extracted 

Real  Part 

Imaginary  Part 

100 

-1.5734  ± j 3.2773 

11.45 

3.64 

-2.1637  ± j 3.8707 

14.97 

15.06 

200 

-1.9994  ± j 3.1408 

0.017 

0.078 

-2.9970  ± j 3.7016 

0.053 

0.192 
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of  two  pole  pairs  than  for  only  one  pole  pair.  The  accuracy,  however 
for  200  samples  is  still  very  satisfactory. 
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Example  3 

As  a final  example  consider  the  transient  response  made  up  of  three 
pole  pairs  such  that 

R(t)  ■ e sin  itt  + e sin  1.5  irt  + e sin  2irt  . (7.13) 

This  signal  is  similar  to  those  of  Examples  1 and  2 in  that  it  decays 
to  very  near  zero  at  about  the  fiftieth  time  step.  The  results  of 
applying  the  Pad4  approximation  to  this  response  are  tabulated  in 
Table  7.3.  Note  that  satisfactory  results  are  never  really  obtained 
for  the  imaginary  parts  of  the  second  and  third  pole  pairs,  because 
the  approximations  used  in  the  procedure  cause  the  Pad4  method  to  have 
difficulty  in  extracting  more  than  two  sets  of  poles.  This,  of  course, 
is  a very  limiting  factor  since  in  almost  all  cases  one  would  desire  to 
obtain  at  least  six  pole  pairs. 

7 . 3 Conclusions  Regarding  the  Pade  Approximation 

It  has  been  shown  that  the  Pad4  approximation  performed  satisfactorily 
when  only  one  or  two  pole  pairs  were  desired  but  failed  to  give  satis- 
factory poles  when  three  pole  pairs  were  sought.  The  method  was  also 
applied  to  signals  with  more  than  three  pole  pairs  and  which  did  not  damp 
out  so  quickly,  but  again  the  results  were  unsatisfactory.  For  the  cases 
of  one  and  two  pole  pairs  it  was  necessary  to  Include  essentially  all  of 
the  waveform  until  it  damped  to  zero.  This  in  itself  is  a limiting 
factor.  Thus,  the  conclusion  would  then  be  that  the  Pad4  approximation,  as 
it  was  used  here,  is  not  a satisfactory  method  for  extracting  the  true 
poles  from  a set  of  transient  response  data. 
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TABLE  7.3 

THE  RESULTING  POLES  USING  THE  PACE  APPROXIMATION  ON  THE  TRANSIENT 
RESPONSE  OF  EXAMPLE  3:  R(t)  - e~2t  sin  irt  + e"3t  sin  1.5  irt  + e“4t  sin  2 irt 


Poles  Extracted 

Percent 

Error 

Samples  Used 

Real  Part 

Imaginary  Part 

-0.9671  ± j 0.2879 

27.73 

76.62 

100 

-2.0313  ± j 3.1523 

17.34 

27.93 

-3.5884  ± j 5.3760 

5.53 

12.18 

-1.9991  ± j 3.1404 

0.023 

0.032 

250 

-2.9915  ± j 4.5577 

0.152 

2.768 

-4.0244  ± j 6.0014 

0.328 

3.783 

-1.9991  ± j 3.1404 

0.023 

0.032 

400 

-2.9915  ± j 4.5577 

0.153 

2.768 

-4.0244  ± j 6.0014 

0.327 

3.782 

8.  CONCLUSIONS 


The  numerical  methods  described  in  this  thesis  provide  a technique 
by  which  the  true  complex  natural  resonances  of  a system  may  be  extracted 
from  a set  of  discrete  transient  response  data  of  that  system.  The 
numerical  procedure,  known  as  Prony's  algorithm,  is  applicable  to  systems 
possessing  multiple  poles  as  well  as  simple  poles.  There  are  two  ef- 
ficient and  systematic  methods  by  which  it  is  possible  to  determine  the 
number  of  poles  inherent  in  the  transient  response.  In  employing  both  of 
these  schemes  the  poles  are  obtained  as  a by-product  of  the  calculations. 
The  Householder  orthogonalization  method,  which  is  the  most  systematic 
of  the  two  procedures  for  determining  the  number  of  poles,  breaks  down 
if  any  significant  noise  is  added.  The  eigenvalue  procedure  on  the  other 
hand  uses  the  known  standard  deviation  of  the  noise  to  aid  in  determining 
the  number  of  poles  but  is  not  quite  as  systematic  as  the  orthogonalization 
procedure. 

The  problem  of  noise  in  the  transient  response  is  a critical  point. 

As  mentioned  above,  the  standard  deviation  of  the  noise  is  used  to  aid 
in  determining  the  number  of  poles.  However,  Prony's  method  will  not 
extract  the  true  poles  if  the  signal-to-noise  ratio  is  below  a certain 
level.  This  level  appears  to  be  in  the  10  to  20  dB  range  for  the  in- 
dividual signal  components  comprising  the  entire  signal.  The  signal-to- 
noise  level  for  the  total  signal  needs  to  be  as  good  as  30  dB.  This 
fact  causes  Prony's  method  to  be  limited  to  use  with  extremely  clean 
data  systems.  However,  if  several  sets  of  noisy  transient  responses  are 
measured  and  averaged  together,  the  standard  deviation  of  the  noise 
decreases  as  one  over  the  square  root  of  the  number  of  trials  run. 
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This  approach  could  be  applied  to  laboratory  test  systems  where  a test 
can  be  repeated  as  many  times  as  desired.  More  work  should  be  done  with 
actual  noisy  experimental  data  where  the  data  are  measured  specifically 
for  the  application  of  Prony's  method.  The  results  would  then  allow 
for  firm  conclusions  to  be  made  about  the  influence  of  experimental 
noise  on  Prony's  method. 

There  are  several  applications  of  Prony's  method  including  system 
analysis,  radar  target  recognition,  the  study  of  spectral  characteristics, 
and  data  reduction  and  extrapolation.  Prony's  method  is  an  invaluable 
tool  in  data  reduction  and  in  the  determination  of  the  spectral  charac- 
teristics of  a system.  The  application  of  the  method  to  radar  target 
recognition  is  also  very  exciting  but  seems  to  be  limited  by  the  noise  level 
problem.  More  work  should  be  done  in  applying  Prony's  method  to  the 
areas  mentioned  as  well  as  to  new  areas. 

The  Pad£  approximation  is  an  alternative  to  Prony's  method.  The 
study  of  this  method  shows  that  the  Padd  approximation  is  not  applicable 
to  systems  containing  more  than  one  or  two  pole  pairs.  This  restricts  its 
usefulness  for  SEM  related  problems.  No  other  alternatives  to  Prony's 
method  were  found. 

The  whole  field  of  SEM  is  new  and  is  growing  so  fast  that  the  outlook 
for  use  of  the  techniques  presented  and  developed  in  this  thesis  is  very 
promising.  It  is  hoped  that  these  techniques  will  help  simplify  the 
study  and  the  analysis  of  the  more  complicated  electromagnetic  structures. 
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